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Executive  Summary 


Classification  and  Discrimination  of  Sources  with  Time-Varying 
Frequency  and  Spatial  Spectra 

ONR,  Grant  no.  N000 14-98- 1-0 176 
Moeness  Amin  (PI) 

This  report  presents  the  results  of  the  research  performed  under  ONR  grant  number 
NOOO 14-98- 1-0 176  over  the  period  of  October  1st,  01  to  September  30th,  02.  The  research  team 
working  on  this  project  consists  of:  Prof.  Moeness  Amin  (PI),  Dr.  Gordon  Frazer  (DSTO, 
Australia),  Dr.  Yimin  Zhang  (Postdoctral  Fellow),  Dr.  Fauzia  Ahmad  (Postdoctoral  Fellow),  Mr. 
Behzad  Mohammadi  (Graduate  Student).  We  have  also  collaborated  with  Dr.  Alex  Gershman 
from  McMaster  University,  Canada,  Prof.  H.  Ge  from  NJIT,  Dr.  N.  Ma  from  DSO,  Singapore, 
and  Prof.  A.  Zoubir  from  Curtin  University,  Australia.  The  research  over  this  fiscal  year  has 
produced  6  journal  papers  and  8  conference  papers.  Copies  of  the  principle  publications  are 
included  in  the  report. 

The  research  efforts  overt  the  past  year  have  focused  on  improved  characterization, 
separations,  suppression,  discrimination,  and  localization  of  nonstationary  and  cyclostationary 
sources  using  multi-sensor  array  receivers.  The  major  contributions  over  the  2001-2002  fiscal 
year  are:  1)  Characterization  of  near-field  scattering  using  quadratic  sensor-angle  distributions,  2) 
Improved  blind  separations  of  nonstationary  sources  based  on  spatial  time-frequency 
distributions,  3)  Introducing  a  unified  representation  of  nonstationary  and  cyclostationary  signals, 
4)  Devising  aperture  synthesis  for  a  through-the-wall  imaging  system,  5)  Nonstationary 
interference  suppression  in  direct  sequence  spread  spectrum  communications  using  space-time 
oblique  projection  techniques,  6)  Formulating  the  time-frequency  ESPRIT  for  direction-of-arrival 
estimation  of  chirp  signals,  7)  Automatic  classification  of  auto-  and  cross-terms  of  time- 
frequency  distributions  in  antenna  arrays,  8)  Proposing  a  new  approach  to  jammer  suppression  for 
digital  communications. 


1.  Characterization  of  Near-Field  Scattering  Using  Quadratic  Sensor-Angle 
Distributions 

We  have  introduced  the  quadratic  sensor  angle  distribution  (SAD)  for  near-field  source 
characterization.  The  SAD  is  a  joint-variable  distribution  and  a  dual  in  sensor  number  and  angle 
to  Cohen's  class  of  time-frequency  distributions.  It  provides  the  power  at  every  angle  for  each 
sensor  in  the  array.  In  this  distribution,  near-field  sources  have  different  angle  for  each  sensor. 
We  use  a  known  test  source  to  illuminate  the  local  scatterer  distribution  we  wish  to  characterize. 
The  high-power  test  source  can  be  removed  via  orthogonal  projection  so  as  to  reveal  the  less 
powerful  local  scatter.  It  is  shown  that  the  eigen-decomposition  of  the  quadratic  representation  of 
SAD  lends  itself  to  source  representation  via  multiple  subarray  beamforming.  The  SAD  can  be 
used  to  clearly  identify  scatterers  on  the  array  axis  both  within  and  beyond  the  array  extent. 
Distinction  between  far-field  spatial  spread  source  and  near-filed  point  source  can  also  be  easily 
achieved. 
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2.  Improved  Blind  Separations  of  Nonstationary  Sources  Based  on  Spatial  Time- 
Frequency  Distributions 

Blind  source  separation  based  on  spatial  time-frequency  distributions  (STFDs)  provides 
improved  performance  over  blind  source  separation  methods  based  on  second-order  statistics, 
when  dealing  with  nonstationary  signals  that  are  localizable  in  the  time-frequency  (t-f)  domain. 
In  the  existing  STFD  methods,  the  covariance  matrix  is  first  used  to  whiten  the  data  vector,  then 
the  mixing  matrix  and  subsequently  the  source  waveforms  are  estimated  using  STFD  matrices 
constructed  from  the  source  t-f  autoterm  points.  We  have  improved  the  STFD-based  source 
separation  method  by  performing  both  whitening  and  estimation  steps  using  the  source  t-f 
signatures.  This  modification  provides  robust  performance  to  noise,  and  allows  reduction  of  the 
number  of  sources  considered  for  separation. 

3.  A  Unified  Representation  of  Nonstationary  and  Cyclostationary  Signals 

The  cyclic  auto-correlation  function,  commonly  used  for  cyclostationary  signals,  and  the 
ambiguity  function,  typically  employed  for  analysis  and  recovery  of  nonstationary  signals,  such 
as  FM,  have  the  same  formulation.  However,  nonstationary  and  cyclostationary  signals  have 
distinct  localization  properties  in  the  time-lag  frequency-lag  domain.  Therefore,  nonstationary 
and  cyclostationary  signals  can  be  represented  and  processed  within  the  same  framework  for 
many  applications,  with  the  distinct  signatures  allowing  effective  source  discriminations.  An 
example  in  array  processing  is  given  where  nonstationary  and  cyclostationary  signals  are 
separated  following  simple  spatial  signature  estimation  exploiting  the  aforementioned  properties. 

4.  Aperture  Synthesis  for  a  Through-the-Wall  Imaging  System 

An  aperture  synthesis  scheme  using  subarrays  that  is  based  on  the  coarray  concept  is 
proposed  for  through-the-wall  microwave  imaging  applications.  Simulation  results  depicting  the 
effectiveness  of  the  proposed  synthetic  aperture  technique  for  a  TWI  system  are  presented.  The 
effects  of  incorrect  estimates  of  wall  parameters  and  errors  in  array  element  placement  on  this 
technique  are  also  investigated. 

5.  Nonstationary  Interference  Suppression  In  DS/SS  Communications  Using  Space- 
Time  Oblique  Projection  Techniques 

Orthogonal  space-time  subspace  using  a  multi-antenna  array  has  been  recently 
proposed  as  an  effective  means  for  nonstationary  interference  suppression  in  direct- 
sequence  spread-spectrum  (DS/SS)  communications.  Interference  mitigation  techniques 
may  compromise  the  desired  signal  in  the  process  of  removing  the  interferers.  This  effect 
is  known  as  self-noise.  The  orthogonal  subspace  projection  introduces  self-noise  and 
signal  distortions  due  to  potentially  high  correlation  between  the  spatio-temporal 
signature  of  the  DS/SS  signal  and  that  of  the  interferers.  We  propose  to  use  the  space- 
time  oblique  projection,  instead  of  the  orthogonal  subspace  projection.  The  oblique 
projection,  while  completely  suppresses  interferers,  does  not  distort  the  desired  signal 
and,  therefore,  no  self-noise  is  generated.  The  performance  of  a  receiver  with  space-time 
oblique  projection,  along  with  the  sensitivity  to  the  errors  in  spatio-temporal  signatures  of 
the  signals  and  interferers,  is  analyzed  and  compared  with  the  receiver  performance 
employing  the  orthogonal  subspace  projection  method. 
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6.  Time-Frequency  ESPRIT  for  Direction-of-Arrival  Estimation  of  Chirp  Signals 

We  have  developed  an  ESPRIT-type  algorithm  for  estimating  the  Directions-Of-Arrival 
(DOA's)  of  multiple  chirp  signals  using  Spatial  Time-Frequency  Distributions  (STFD's).  An 
averaged  STFD  matrix  (or  multiple  averaged  STFD  matrices)  are  used  instead  of  the  covariance 
matrix  to  estimate  the  signal  and  noise  subspaces.  The  proposed  algorithm  is  shown  to  provide  a 
significant  performance  improvement  over  the  traditional  ESPRIT  algorithm  for  FM  sources, 
specifically  in  situations  with  closely  spaced  sources  and  low  Signal-to-Noise  Ratios  (SNR's). 
Simulation  results  are  provided  to  illustrate  the  performance  of  the  proposed  approach  in 
scenarios  with  multiple  narrowband  chirp  sources. 

7.  Automatic  Classification  of  Auto-and  Cross-Terms  of  Time-Frequency 
Distributions  in  Antenna  Arrays 

The  problem  of  selecting  auto-  and  cross-terms  of  time-frequency  distributions  (TFDs)  of 
nonstationary  signals  impinging  on  a  multi-antenna  receiver  is  considered.  A  detection  approach 
is  introduced  which  allows  performance  measurement  and  comparison  of  various  schemes  via 
receiver  operating  characteristics.  Array  averaging  and  array  differencing  techniques  are  both 
employed  to  form  a  basis  for  time-frequency  point  selection.  The  proposed  classification  method 
is  evaluated  against  the  bootsrap-based  method.  It  is  shown  that  the  former  offers  improved 
performance  and  simplified  implementations. 

8.  A  New  Approach  to  FM  Jammer  Suppression  for  Digital  Communications 

We  have  considered  the  problem  of  FM  jammer  suppression  in  digital  communications. 
It  is  pointed  out  that  the  cyclic  auto-correlation  function,  commonly  used  for  cyclostationary 
signals,  and  the  ambiguity  function,  typically  employed  for  analysis  and  recovery  of 
nonstationary  signals,  such  as  FM,  have  the  same  formulation.  However,  nonstationary  and 
cyclostationary  signals  have  distinct  localization  properties  in  the  time-lag  frequency-lag  domain, 
and  this  property  can  be  effectively  used  for  jammer  suppression.  Utilizing  the  spread  of  an  FM 
jammer  signature  beyond  that  of  the  communication  signal,  one  can  select  the  jammer-only  time- 
lag  frequency-lag  points  for  effective  jammer  spatial  signature  estimation.  Suppression  of  the 
jammer  signal  is  then  proceeded  by  projecting  the  received  signal  to  the  jammers'  orthogonal 
spatial  subspace.  Simulation  examples  for  the  recovery  of  BPSK  signals  in  the  presence  of  strong 
and  moderate  FM  jammers  are  presented. 
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Characterization  of  Near-Field  Scattering  Using 
Quadratic  Sensor- Angle  Distributions 


Gordon  J.  Frazer  and  Moeness  G.  Amin 


Abstract 

We  introduce  the  quadratic  sensor  angle  distribution  (SAD)  for  near-field  source  characterization. 
The  SAD  is  a  joint-variable  distribution  and  a  dual  in  sensor  number  and  angle  to  Cohen’s  class  of 
time-frequency  distributions.  It  provides  the  power  at  every  angle  for  each  sensor  in  the  array.  In 
this  distribution,  near-field  sources  have  different  angle  for  each  sensor.  We  use  a  known  test  source 
to  illuminate  the  local  scatterer  distribution  we  wish  to  characterize.  The  high-power  test  source  can 
be  removed  via  orthogonal  projection  so  as  to  reveal  the  less  powerful  local  scatter.  It  is  shown  that 
the  eigen-decomposition  of  the  quadratic  representation  of  SAD  lends  itself  to  source  representation  via 
multiple  subarray  beamforming.  The  SAD  can  be  used  to  clearly  identify  scatterers  on  the  array  axis 
both  within  and  beyond  the  array  extent.  Distinction  between  far-field  spatial  spread  source  and  near-filed 
point  source  can  also  be  easily  achieved. 
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I.  Introduction 


In  this  paper,  we  propose  a  new  method  for  characterizing  the  near-field  electromagnetic 
scatter  local  to  a  linear  equi-spaced  array  of  electromagnetic  sensors.  The  particular  application 
is  from  the  area  of  high  frequency  surface-wave  radar,  although  the  proposed  method  may  be 
also  applied  in  other  sensor  array  applications,  provided  the  array  is  linear  and  equi-spaced  and 
the  local  scatterers  are  in  the  near-field. 

Surface-wave  radar  is  an  emerging  coastal  and  exclusive  economic  zone  surveillance  technol¬ 
ogy.  Surface-wave  radar  systems  used  in  this  role  are  located  at  a  land-sea  boundary  and  use 
surface-wave  propagation  and  the  conductivity  of  sea  water  to  detect  and  track  targets  across 
water  beyond  the  line-of-sight  radar  horizon.  Detection  and  tracking  of  small  surface  vessels  can 
be  achieved  at  ranges  in  excess  of  200km  where  the  optical  horizon  may  be  no  more  than  20km. 
These  radars  operate  in  the  congested  lower  HF  (approximately  3-10MHz)  section  of  the  elec¬ 
tromagnetic  spectrum.  The  systems  with  which  we  have  experience  use  floodlight  transmission 
and  an  array  of  receiver  sensors.  A  mix  of  classical  digital  beamforming,  space-time  adaptive 
processing,  and  high  resolution  angle  algorithms  are  used  to  determine  target  direction. 

A  typical  surface- wave  radar  receiving  array  may  consist  of  between  8  and  64  sensors  and  can 
be  500m  or  1km  in  total  length.  It  is  typically  sited  on  a  coastal  beach  which  may  or  may  not 
provide  a  uniform  transition  from  land  to  sea.  For  example,  the  coast  may  in  fact  be  a  bay  in 
which  case  the  land  sea  boundaries  beyond  either  end  of  the  array  may  cause  near-field  scattering 
and  distort  the  wave-front  arriving  at  the  array.  There  may  be  other  locally  sited  structures,  such 
as  buildings  and  fences,  which  can  be  the  source  of  local  scatter  (consider  that  the  wavelength 
of  the  radar  signal  is  between  30-100m).  This  makes  achieving  very  low  sidelobe  spatial  beams 
with  a  classical  beamformer  a  difficult  problem  and  can  render  the  receiver  system  vulnerable  to 
interference  through  beam  sidelobes  (possibly  via  sky  wave  propagation). 

The  near-field  scatter  produced  by  these  mechanisms  are  correlated  with  the  desired  direct 
far-field  radar  return  from  targets  (and  clutter).  This  scatter  is  typically  approximately  20- 
40dB  weaker  than  the  direct  signal.  Without  mitigation  it  is  possible  to  achieve  classical  beam 
sidelobes  of  30-35dB,  however  in  general  the  remaining  components  of  the  receiving  system  can 
sustain  substantially  higher  performance  [1]. 

A  method  is  required  that  can  mitigate  the  effect  of  local  scatter  so  that  the  radar  system  can 
realize  the  inherent  sidelobe  capability  as  set  by  the  radar  equipment  [1].  An  important  part  of 
this  strategy  is  to  first  characterize  the  local  scatter  distribution.  A  means  of  performing  this 
characterization  is  the  subject  of  the  work  we  present  in  this  paper. 

Breed  and  Posch  introduced  the  spatial  Wigner  distribution  (SWD)  as  a  tool  for  determining 
the  range  and  angle  of  a  near  field  source  [2].  They  exploited  the  property  that  the  phase  front 
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of  a  wave  emanating  from  a  source  in  the  near  field  of  an  array  has  an  approximately  quadratic 
phase  law,  or  equivalently  an  approximately  linear  spatial  frequency  law.  The  true  propagating 
wave  phase  front  is  in  fact  spherical  and  is  only  approximately  quadratic  for  near  field  sources 
some  distance  from  the  array,  and  therefore,  the  method  in  [2]  breaks  down  for  sources  close  to 
the  array. 

A  comprehensive  treatment  of  SWD  was  given  by  Swindlehurst  and  Kailath  in  [3].  They 
included  an  examination  of  the  applicability  of  the  quadratic-spherical  approximation  and  use 
a  parametric  high  resolution  technique  to  determine  the  linear  frequency  law  parameters  (and 
hence  the  near-field  source  position). 

More  recently  a  related  but  altogether  different  spatial  time-frequency  distribution  (STFD) 
was  introduced  in  [4],  [5].  This  is  a  true  spatial  TFD  in  that  individual  entries  in  the  array 
spatial  covariance  matrix  are  replaced  by  time-frequency  representations  of  the  energy  comprising 
these  entries  (both  auto  and  cross)  throughout  the  interval  T.  This  approach  has  been  applied 
successfully  to  two  challenging  array  processing  problems;  blind  signal  separation  and  angle 
estimation. 

There  is  a  substantial  body  of  literature  concerned  with  processing  spatial  signals  received  by 
an  array  of  sensors  from  sources  in,  the  near-field  of  the  array,  e.g.  [6].  It  is  mostly  concerned 
with  techniques  for  estimating  the  angle  and  range  of  the  source.  Both  a  subspace  method  and 
a  maximum  likelihood  algorithm  are  presented  in  [7] . 

Several  authors  have  proposed  methods  for  determining  the  angle  of  distributed  sources  lo¬ 
cated  in  the  far-field  of  an  array  [8],  [9].  These  techniques  address  the  effect  of  scatter  local  to  a 
transmitter  in  the  far-field  and  not  for  scatter  that  is  sufficiently  local  to  the  receiving  system  to 
be  in  the  near  field  of  the  array.  In  this  paper,  we  propose  a  generalization  of  the  spatial  Wigner 
distribution  introduced  in  [2],  combined  with  orthogonal  projection  techniques  for  the  character¬ 
ization  of  the  structure  of  local  electromagnetic  scattering  induced  by  the  near  environment  and 
mutual  interaction  between  sensors  in  the  array.  We  use  a  known  test  source  to  illuminate  the 
local  scatterer  distribution  we  wish  to  characterize.  The  high-power  test  source  can  be  removed 
via  orthogonal  projection  so  as  to  reveal  the  less  powerful  local  scatter.  It  is  shown  that  the 
eigen-decomposition  of  the  quadratic  representation  of  SAD  lends  itself  to  source  representation 
via  multiple  subarray  beamforming.  The  SAD  can  be  used  to  clearly  identify  scatterers  on  the 
array  axis  both  within  and  beyond  the  array  extent.  Distinction  between  far-field  spatial  spread 
source  and  near-filed  point  source  can  also  be  easily  achieved. 

This  paper  is  organised  as  follows.  In  section  2,  a  signal  model  for  far-field  and  near-field 
sources  are  presented.  Both  point  and  spatial  spread  sources  are  considered.  The  sensor-angle 
distribution  is  introduced  in  section  3.  The  source  range  and  angle  expressions  viewed  by  each 
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sensor  are  presented  in  section  4  using  the  geometrical  relation  between  the  source  of  interest 
and  the  array.  Exact  and  least-squares  maximum  estimates  of  the  source  range  and  angle, 
along  with  Cramer-Rao  bound,  measured  from  the  center  of  the  array  using  the  respective  SAD 
estimates  at  each  sensor  are  also  derived  in  section  4.  Alias-free  implementations  and  subarray 
beamforming  interpretation  of  the  SAD  are  delineated  in  Section  5.  Section  6  includes  simulations 
demonstrating  the  offerings  of  the  SAD  in  near-field  source  characterization. 

II.  Signal  Model 

A.  Model 

Our  proposed  measurement  technique  requires  one  cooperative  source,  S,  with  complex  enve¬ 
lope  s®,  in  the  far-field  of  the  array  at  known  angle  6s  and  where  k  denotes  time  index.  Consider 
a  linear  equi-spaced  array  of  M  sensors,  where  sensor  position  errors  are  negligible  and  the  gain 
and  phase  of  all  sensors  are  accurately  matched.  Steering  vectors  for  the  far-field  a(0)  and  near¬ 
field  a(0,  r)  take  on  the  standard  form  with  6  being  the  angle.  For  the  near-field  steering  vector, 
r  denotes  range  [10]. 

Assume  that  the  conditions  on  the  test  source  and  sensor  array  are  such  that  the  following 
signal  model  is  appropriate 

zk  =  ASk  +  qk  +  nk  (1) 

In  this  model,  zk  is  a  vector  of  dimension  M,  representing  the  kth  snapshot  of  sensor  outputs.  The 
vector  qk  consists  of  additive  spatial  and  temporal  coloured  noise  produced  in  the  environment, 
and  nk  represents  additive  white  noise  modeling  the  internal  noise  of  the  array  of  sensors  receiving 
system. 

The  matrix  A  can  take  on  one  of  two  forms,  depending  on  whether  the  local  scatterer  is  best 
modeled  as  a  collection  of  P  discrete  point  scatterers,  or  as  a  single  spatially  distributed  scatterer. 
For  the  case  of  P  discrete  scatterers 

A  =  [a(0s),a(0i,r1)>...  ,a(0i, ,a(0P,rP)]  (2) 

In  the  above  equation,  the  near-field  scatterer  i  =  1, . . .  ,  P  is  characterized  by  the  angle  0j  and 
range  r,-.  In  the  case  of  a  distributed  scatterer  contained  in  the  near  field  azimuth  and  range  set 
Q,  with  reflectivity,  i?n(0,r)  and  with  respect  to  a  reference  location  (0o ,  ro ) 

n  /nRp(0;  —  0p,r'  —  r0)  ■  a(0'  —  0p,r'  —  r0)  d0'dr' 

a  1  0,r°'  ||/nRn(0'-0o,r'-ro)a(0'-0o,r'-ro)  d0'dr'||  1  ’ 

and  hence 

A=[a(0s),afi(0o,ro)]  (4) 
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Likewise,  the  signal  vector  Sk  can  be  constructed  in  two  ways,  depending  on  whether  the 
near-field  scatterers  are  best  modeled  as  discrete  or  distributed.  For  the  discrete  case, 

Sk=  [sk,Sk,-.:  ,Sk]  (5) 

where  the  test  source  complex  amplitude  is  given  by  and  the  complex  amplitude  of  the  ith  of 
P  discrete  scatterers  is  denoted  as  s\.  In  the  distributed  scatterer  case, 

Sk=[4,Sk]  (6) 

The  and  s\  and  may  be  uncorrelated  or  correlated  for  each  case  respectively. 

The  ith  element  of  the  steering  vector  for  the  far-field  source  (the  first  column  of  A)  a (0s) 
takes  on  the  standard  form 


a.i(8)  =  exp  -j 


.2iui  . 


(i-1)- 


M-  1 


whereas,  in  the  case  of  any  near-field  source  or  scatterer,  the  ith  element  of  the  steering  vector 
a(0s,  rs),  is  given  by 

ai(rs ,  0S)  -  —  •  exp  ( -j ^  •  rS)A  (8) 


ai(rs,Vs)  -  —  •  exp  •  rs,i J  (8) 

assuming  a  normalized  and  equal  gain  for  each  sensor. 

The  spatial  environment  is  characterized  by  the  spatial  covariance  matrix  Rk,m>  where 

Rk,m  =  AS^mA  +Qk,m"t~cr  I  (9) 

with  Rk,m  =  E[zkZ^],  the  source  covariance  Sk,m  —  E[skS**],  and  the  noise  covariance  Qk,m  — 
Assume  that  Sk,  qk  ^.nd  ilk  3-re  uncorrelated,  that  S^m  —  S^k—m  ^nd  Qk,m  —  Q^k— mi 
and  that  S5k  =  S<Jk'  and  Q^k  =  Q^k'  for  k  different  from  kl .  In  this  case  S,  Q  and  hence  R 
are  temporally  white  and  stationary,  and  as  such  we  can  remove  the  dependency  on  k  and  m. 
These  assumptions  are  not  in  fact  strictly  required  but  we  include  them  in  order  to  reasonably 
bound  the  class  of  signals  we  consider  in  the  subsequent  discussions.  The  variance  of  the  receiver 
internal  noise  is  a2.  Individual  elements  of  matrix  S  are  denoted  as  pij. 

We  ensure  that  the  cooperative  test  source  has  sufficient  signal  to  noise  ratio  (generally  greater 
than  50dB)  to  perform  our  measurement  by  requiring  that 

Psnr  =  W  +  tr[Q])  >>  1  (10) 

It  is  also  expected  that  the  direct  far-field  source  power  will  be  substantially  greater  than  the 
total  near-field  power  (by  20-40dB),  i.e, 


Psni  = 
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III.  Sensor-Angle  Distributions 

The  proposed  method  uses  an  extension  to  the  spatial  Wigner  distribution  originally  introduced 
by  Breed  and  Posch  [2],  To  avoid  confusion,  it  has  been  necessary  to  change  the  name  to  reflect 
the  generalization  to  all  members  of  Cohen’s  class  of  quadratic  distributions  [11].  While  the 
title  “spatial  Wigner  distribution”  is  informative,  retaining  the  name  “spatial  time-frequency 
distribution”  for  the  remaining  members  of  Cohen’s  class  applied  to  spatial  signals  does  not 
correctly  describe  the  distribution  we  are  interested  in.  Therefore,  in  this  work,  we  have  renamed 
the  class  of  quadratic  distributions  applied  to  spatial  signals  to  be  sensor-angle  distributions 
(SAD).  The  corresponding  spectra  are  called  sensor-angle  spectra  (SAS). 

The  Cohen’s  class  SAD  for  the  kth  snapshot  is  a  distribution  of  the  angle  of  sources  impinging 
on  the  array  at  each  sensor , 

oo  oo 

Tk(i,0;zk)=]T  Y  <Km>1)Mi  +  m  +  1)4(i  +  m-1)e_j4^1  (12) 

1=— oo  m=— oo 

where  i  and  6  are  the  sensor  index  and  angle  respectively.  The  kernel  0(m,  l)  characterizes  the 
distribution  and  is  a  function  of  sensor  position,  m,  and  sensor  lag,  l.  All  the  standard  data- 
independent  or  data-dependent  kernel  designs  applied  in  the  time-frequency  literature  may  be 
used  with  the  SAD  [11]. 

The  sensor-angle  spectrum  (SAS)  is  the  power  (not  energy  or  energy  density)  distribution  of 
the  sources  impinging  on  the  array.  The  SAS  is  given  by 

Ts(i,0;zk)=E[Tk(i,0;zk)]  (13) 

where  an  estimate  for  temporally  stationary  sources  is  given  by 

1  N_1 

fS(i,0;zk)  =  -^Tk(i,0;zk)  (14) 

k=0 

for  N  snapshots. 

The  objective  is  to  the  SAD  (or  SAS)  to  characterise  the  near-field  scatterers  of  a  far-field 
source  signal.  We  expect  the  test  signal  to  be  substantially  more  powerful  than  the  local  scatter 
we  wish  to  characterize  (see  (11)).  With  the  knowledge  of  the  far-field  source  angle-of-arrival, 
a  spatial  filter  can  be  designed  to  remove  its  dominance,  allowing  a  clear  depiction  of  the  near 
field  source  in  the  sensor-angle  (s-a)  domain.  Alternatively,  a  simple  technique  is  to  project  the 
sensor  data  on  the  orthogonal  subspace  of  the  far-field  source  spatial  signature.  In  (12)  and  (13), 
the  data  snapshot  zk  is  replaced  by  Pzk  where  P  is  the  orthogonal  projection  operator  formed 
from  the  test  source  steering  vector  a (0s)  as 

P  =  I  -  a(0s)[aH(0s)a(0s)]-1aH(0s)  (15) 
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Therefore,  we  compute  the  SAS 


Ts(i,0;P)|,s  (16) 

associated  with  a  test  source  in  the  far-held  of  the  array  at  angle  0s. 

In  some  applications,  a  single  test  angle  will  provide  sufficient  characterization  using  (16) 
while  in  other  applications,  two  or  many  test  angles  will  be  required,  in  which  case  the  full 
characterization  is  given  by 

Ts(i,0;P(0s))  (17) 

as  0s  is  scanned  over  the  required  domain  of  angles  for  the  test  source. 


IV.  Range-Angle  Estimation  from  SAD 

A .  Geometry 
A.l  Assumptions 

Consider  a  linear  equi-spaced  array  of  M  sensors  placed  on  a  flat  plane  in  a  two  dimensional 
surface.  Assume  that  sensor  position  errors  are  negligible  and  the  gain  and  phase  of  all  sensors 
and  corresponding  data  acquisition  equipment  are  accurately  matched.  Assume  that  the  array 
is  narrowband,  i.e.,  the  reciprocal  of  the  bandwidth  of  any  signals  received  is  large  compared 
with  the  propagation  delay  across  the  array.  The  wavelength  of  all  sources  received  is  A.  Let  the 
origin  of  a  coordinate  system  O  be  at  the  mid-point  of  the  array,  with  the  sensors  individually 
spaced  by  d  regularly  along  the  x-axis  and  indexed  i  =  ,M  from  left  to  right.  We  assume 
that  d  <  |.  Boresight  is  along  the  y-axis. 

A  source  is  placed  in  the  near-field  (i.e.  a  circular  wavefront  impinges  on  the  array)  at  location 
rs  meters  from  the  origin  and  0S  degrees  from  boresight.  We  define  that  angles  are  to  be  measured 
clockwise  from  array  boresight  (the  y-axis).  For  M  odd,  there  is  a  sensor  at  the  origin,  whereas 
for  M  even,  the  origin  is  midpoint  between  two  sensors.  The  array  geometry  and  the  notations 
are  shown  in  figure  1  for  the  case  of  M  =  8. 


A. 2  Sensor  angle  for  a  source  at  (rs,  6S) 

The  physical  position  of  some  point  along  the  axis  of  the  array,  with  respect  to  an  origin  O  at 
the  midpoint,  is  denoted  x.  The  position  of  the  ith  sensor  in  the  array  can  be  written 


x  =  d  • 


(*-!)“ 


for  i  G  , M} 


The  distance  from  position  a:  to  a  source  at  (rs,0s)  is 

rSyX  =  \fx2  -  2  ■  rs  ■  x  ■  smQs  +  r2s 


(18) 


(19) 
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with  the  angle  given  by 


—  cos 


X2  +  rlx  -  rj 

2  •  x  ■  rS)I 


7T 

2 


(20) 


Substituting  equation  (19)  into  equation  (20)  eliminates  the  dependency  on  rs>x  so  that  0SjX, 
which  is  the  angle  from  position  x  to  the  given  source  at  (rs,6s),  may  now  be  written  directly  in 
terms  of  the  source  location 


6s,x  =  cos 


x  —  rs  •  sin  6S 

\fx*  —  2  •  rs  •  x  •  sin  0S  +  rf 


(21) 


A  specific  case  of  interest  is  the  angle  at  each  sensor,  which  may  be  expressed  as  a  function  of 
the  source  location  and  the  sensor  number  i  €  ,  M}. 


0s  i  =  0(z,rs,0s)  for  a  given  d  and  M 


(22) 


Substituting  equation  (18)  into  equation  (21)  gives  a  direct  expression  for  the  angle  at  the  ith 
sensor  to  a  source  at  location  (rs, #,) 


@(i,rs,6s)  =  cos" 


[d  ■  [(i  -  1)  -  ^j1]]  ~  rs  ■  sin 6S 


.\f[d  ■  [(*  -  1)  -  ^l]2  -  2  •  rs  •  [d  ■  [(*  -  1)  -  ^]]  •  sin0s  +  rf 


7 r 
2 


(23) 


A.3  Source  location  from  sensor  angle 

Given  knowledge  of  the  sensor-angle  at  any  two  or  more  sensors,  it  is  possible  to  determine  the 
range  and  bearing  (rs,05)  of  a  source.  This  is,  however,  subject  to  identifiability  requirements 
that  each  sensor  has  a  different  angle  to  the  source,  which  is  in  fact  a  requirement  that  the  source 
be  in  the  near-field  of  the  array.  Given  6Sji  and  6sj  with  i  /  j,  one  determines  sensor  to  source 
ranges  rSfi  and  rsj  respectively  using 

cos[6Sj] 


rs,i  =  [|*  -  3 1]  ■  d  ■ 


and 


-  IT 


sin[05jj  -  9s,j\ 
cos[g5|j] 


(24) 


(25) 


sin [6Sii  -  es,j] 

This  requires  that  dsj  —  9sj  ^  0.  The  range  and  bearing  with  respect  to  the  origin  can  be 
determined  relative  to  any  of  the  individual  sensors  using  the  individual  sensor  range  and  bearing. 
For  example,  for  the  jth  sensor,  we  use  rSJ  and  9Sij  according  to 


2  2  , 
rs  ~  rs,j 


M  —  1 

2 

—  1 

[— 2 - \j-tf\-d 

—  2  •  rs,j  • 

i — a — b'-in-rf 

sin  0, 


(26) 
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and 


6S  =  cos  1  cos  6sj] 

rs 


(27) 


A.4  Estimating  source  location 

Alternatively,  an  estimate  for  source  location  can  be  determined  given  estimates  of  the  sensor 
angle  0sj  and  equation  (23).  The  least  squares  estimate  for  source  location  is 


=  argmin 

rs  ,0S 


—  ©(*,TS,0S)]2 


and,  under  the  assumption 

[^s,i  ~  Os,i]  “►  -A f  (0,  of) 

then  the  maximum  likelihood  estimate  for  source  location  is 


__  arg  min 

rs,0s 


£[<?s,i-©(i,rs,0s)]2 


(28) 


(29) 


(30) 


B.  Cramer-Rao  Bound  for  Location  from  Sensor  Angle  Measurements 

Estimating  the  angle  to  a  source  at  each  sensor  from  a  spatial  signal  is  analogous  to  estimating 
the  instantaneous  frequency  (IF)  of  a  single  component  time  domain  signal.  Numerous  techniques 
for  estimating  IF  have  been  proposed,  including  methods  based  on  model  fitting  of  unwrapped 
phase  estimates,  and  from  extraction  of  the  peak  value  at  each  time  in  particular  time-frequency 
distributions  [12].  More  recently  a  technique  using  iterative  linear  representations  of  the  signal 
has  been  proposed  based  on  the  cross  polynomial  Wigner  distribution  [13]. 

Performance  bounds  for  IF  estimation  are  given  in  [14]  and  a  comparison  of  several  IF  esti¬ 
mation  techniques  with  respect  to  the  Cramer-Rao  lower  bound  (CRLB)  is  given  in  [13].  These 
results  are  applicable  for  high  signal  to  noise  ratio  cases  and  include  expressions  for  bias  and 
variance  in  the  overall  frequency  function  of  time  estimates  (and  not  at  one  instant  of  time). 
These  results  show  that  the  estimate  will  be  biased  for  IF  laws  other  than  constant  or  linear 
(using  a  second  order  TFD)  or  for  a  polynomial  IF  law  of  order  greater  than  the  order  of  a 
polynomial  TFD. 

Recently,  expressions  for  the  bias  and  variance  of  IF  estimates  based  on  extracting  the  peak 
of  a  TFD  at  each  time  have  been  developed  for  the  case  of  low  signal  to  noise  ratio  (SNR) 
[15].  These  results  explain  the  outlier  behaviour  frequently  observed  when  using  peak  extraction 
algorithms  at  low  SNR.  The  authors  have  partitioned  the  error  behaviour  into  two  mechanisms, 
corresponding  to  high  and  low  SNR,  and  provide  insight  into  when  each  mechanism  applies. 
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Regardless  of  the  particular  IF  estimation  technique,  the  aforementioned  performance  bounds 
governing  IF  estimation  at  each  data  sample  can  be  applied  to  the  problem  of  determining  the 
localization  performance  for  a  near-field  source,  in  range  and  angle,  using  sensor-angle  measure¬ 
ments.  Given  knowledge  of  the  IF  (or  sensor-angle)  estimation  error,  we  now  derive  the  CRLB 
for  the  location  of  the  source  in  range  and  angle. 

The  assumption  that  the  sensor-angle  measurement  estimate  errors  are  Gaussian  and  indepen¬ 
dent  (high  SNR  case)  (equation  (29))  can  be  re-stated  as  the  following  sensor  angle  measurement 
model 


4,i  =  ©(*’  I  rs,  0S)  +  Wi  for  i  E  {1, . . .  ,  M}  (31) 

where  wt  ~  J\f(Q,af)  and  E[w;Wj,]  =  er? <$;„;/  and  where  the  variance  of  Wi  may  be  different  for 
each  sensor  i.  The  assumptions  of  Gaussianity  and  of  independence  of  angle  measurements  at 
different  sensors  are  discussed  in  [16]  in  the  time-frequency  context.  These  assumptions  may  not 
be  satisfied  for  all  members  of  SAD  [15],  however,  they  simplify  the  CRLB  derivations  discussed 
in  the  following  section.  Let 


y  =  [4,1)  •  •  •  Am]13 


be  the  vector  of  sensor  angle  measurements  and  let 


0(rs,0s)  =  [0(l;rs,0s),.. 


,©(M;rs,0s)]T 


(32) 


(33) 


be  the  vector  of  predicted  measurements.  The  likelihood  and  log  likelihood  functions  for  this 
model  are 


M  i  ,  i  . 

p( y ; rs,  es)  =  JJ  - r=  [(9s>i-©(i;rs,0s)]  ) 

i  yj2i raf  \  zai  L  J  / 


(34) 


and 


1  2 

lnp(y ;  rs,  0S)  =  -  ^  ln[V^  ■  a-,}  -  ^  [4,i  -  ©(i ;  rs,  0S)]  (35) 

i  i  * 

Or,  in  vector  form  using  equations  (32)  and  (33)  and  with  E  =  diag(crf, . . .  ,o2M)  (where  diag  is 
the  square  matrix  with  the  specified  entries  along  the  diagonal  and  zero  elsewhere)  we  can  write 
equation  (34)  as 


y  ~  Ar(0(rs,0s),S) 


(36) 
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To  determine  the  CRLB  we  require  the  Fisher  Information  matrix  (FIM)  (for  the  two  parameters 
to  be  estimated  rs  and  6S) 

f  rM02lnp(y  ;rs,6s)i  x?[d2  lnp(y  ^,0,)!  "| 

T  t f\  \  J  - J 


j(r  Q  \  —  L  dri  J  L  0r,08,  J  /oy\ 

S’  S  p,rd2  lnp(y  ;rs,0s)i  pr02  inp(y  ;r5,05)i  ' 

.  "I  d8,dr,  J  "l  gj%  J  _ 

For  the  particular  case  of  a  model  of  the  form  of  equation  (36)  the  FIM  simplifies  to  [17] 

r8e(r.,g,)irv,-irae(r„e,)1  rae(rs,e,)irv-ir^Q(r»,gs)i 
j/r  Q  \  _  1  3r»  J  1  dr,  J  1  Or,  J  1  08,  J  /no\ 

and  S_1  =  diag(^j , . . .  ,  4-).  The  individual  elements  of  I(rs ,6S)  in  equation  (38)  can  be  written 
in  terms  of  the  sensor  elements  as 


x /_  ^  1  ae(i;rs,6>s)  50(i;rs,6»s) 

I(r.,0s)i,i  =  2^3 - - ^ - 


drs 


a  \  1  d©(i;rs,0s)  d©(i;rs,0s) 


t/w  o  \  X^  1  d@{i-rs,es)  de{\-,rs,es) 

^^2,1=  2^-2 - JZ - o- - 


I(rs,0s)2,2  =  E 


1  d©(i;rs,0s)d©(i;rs,0s) 


* — #  a* 
i=X  1 


?  Ms  d6s 


This  requires  that  the  partial  derivatives  of  equation  (23)  with  respect  to  the  location  parameters 
( rs,6s )  be  evaluated 


dQ{i;rs,es) 

drs 


sinfls _ (x—rs  S\n9s)(-2xsm9s+2rs 

-2  rs  x  sin  Os~\~r $ 2  2  (x2-2  rs  x  sin  0s+rs2)3/2 

/l  —  ix~rs  sin  Os)2 
V  1  xz -2  rsx  sin  Os+rs  2 


deji^Qs) 

dOs 


( x-Ta  sin Os)rsx cos 9S 


/x2~2  rsx  sin  Os  +rs2  (s2  -2  rax  sin  0S  +rs 2 ) 3 

t/l  (J~r^  sinfl5)2 

V  1  x2—2ra  xsin 0s+rs2 


with  x  =  d-  [(*  -  1)  -  for  i  G  {1 ,M}  (as  in  equation  (18)).  Let  equation  (37)  be 


re-written  as 


I(rs,  6S)  — 


I(rSA)l,l  I(rs4s)l,2 
I(rs,^s)2,i  I(rs,08)a,2 
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where  the  individual  elements  are  determined  using  equations  (43)  and  (44),  evaluated  at  a 
particular  ( rs,9s )  of  interest,  substituted  into  equations  (39)  through  (42),  where  the  a\  are  the 
variance  of  the  angle  measurement  at  each  sensor.  The  Cramer-Rao  lower  bounds  for  estimates 
of  rs  and  6S  are,  respectively 


var(rs)  > 


I(rs,0s)2,2 


[I(rs,  $s)l,l  ‘  I(rs,0s)2,2  I(rs>  $s)l,2  '  I(rs,  $s)2,l] 


(46) 


and 


var(0s)  > 


I(rs,^s)i,i 


[I(rs,  0s)i,i  •  I(rs,  0S) 2,2  —  I(rs,  0s)i,2  •  I(rs,  0s)2,i] 


(47) 


The  CRLB  is  a  function  of  the  source  location.  For  one  particular  array  configuration  and  sensor- 
angle  estimate  variance  at  each  sensor,  the  variation  of  CRLB  for  source  location  over  a  grid  of 
possible  source  locations  is  shown  in  figures  2  and  3. 


V.  Implementation 

A.  Alias-Free  SAD 

In  practice,  implementing  the  SAD  directly  using  the  previous  definition  has  proven  undesir¬ 
able.  This  can  be  seen  by  rewriting  equation  (12)  as 

Tk(i,0;zk)=JF[KZk(i,l)^(i,l)]  (48) 

where 

KZk(i,l)  =  zk(i  +  \)z*k(i-l)  (49) 

(F  denotes  Fourier  transform  in  the  variable  l  6  and  *  denotes  convolution  in  the  variable  i.) 
The  sensor  position  (index  i )  dependent  lag  (index  l )  sequence,  KZk(i,l),  is  evaluated  at  even 
lag  intervals  only  (i.e.  where  the  lag  spacing  between  sensors  is  even).  This  under-samples  the 
true  position  dependent  lag  sequence  and  causes  aliasing  in  the  resulting  SAD  for  many  source 
locations. 

It  is  possible  to  correct  the  aliasing  problem  by  oversampling  the  spatial  signal  by  two  (i.e. 
space  the  sensors  by  d  =  |)  or  by  interpolating  the  d  =  |  sampled  spatial  signal.  Oversampling 
is  frequently  impractical  as  it  doubles  the  cost  of  the  array.  Interpolation  is  also  undesirable 
because  sensor  arrays  frequently  have  a  limited  number  of  sensors  (M)  and  the  array  “end” 
effects  associated  with  interpolating  the  finite  extent  spatial  signal  corrupt  the  resulting  SAD. 

A  more  satisfactory  approach  is  to  exploit  the  results  of  Jeong  and  Williams  [18].  By  rotating 
the  domain  of  KZy  (i,  l)  and  (f)(i,  l )  by  -45°  it  is  possible  to  construct  values  in  4>(i,  l)  corresponding 
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to  odd  interval  lags  in  KZk{i,l).  Let 

=  zk(ii)zk(i2)  (50) 

and 

=  (51) 

for  sensors  i\  and  *2-  The  rotated  K'Zk  contains  the  |*i  —  h\  G  {•••  ,—2,0,2,...}  lag  values 
corresponding  to  the  even  lags  determined  in  equation  (49).  It  also,  however,  contains  the 
missing  odd  lag  samples  (| ^  G  {. . .  ,  -3,  -1, 1,3,...})  required  to  avoid  aliasing.  Similarly, 
the  kernel  <f>'{i\,i2)  is  evaluated  on  the  rotated  (*1,^2)  domain  and  likewise  provides  support  for 
the  missing  odd  interval  lag  values.  These  values  may  be  interpolated  from  the  even  lag  kernel 
samples,  or,  exactly  computed  from  the  definition  of  the  kernel  sequence.  Note  that  the  rotated 
kernel  (fi'ihih)  has  a  larger  extent  than  prior  to  rotation.  Using  this  approach  it  is  possible  to 
compute  an  approximately  alias-free  SAD  [19]. 

B.  Quadratic  SAD 

Cunningham  and  Williams  [20],  and  Amin  [21],  have  shown  (in  the  context  of  time-frequency 
distributions)  that  it  is  possible  to  express  Tk(i,  0;zk)  in  quadratic  form  as  a  weighted  sum  of 
spatial  beamformers.  Let  the  eigen-decomposition  of  the  rotated  (or  alias-free,  see  equation  (51)) 
kernel  be  given  as 

4>'(h,k)  =  VSVH  (52) 

where  the  columns  vj  (for  j  G  {0, . . .  ,  J  —  1})  of  V  are  the  J  eigenvectors  of  </>'  and  the  diagonal 
entries  of  matrix  S  are  the  corresponding  eigenvalues,  oo> .  .  ■  ,  crj- 1-  Due  to  the  use  of  the  rotated 
kernel  1,^2),  J  =  4L  +  1,  for  a  given  maximum  lag  extent  L  in  the  non-rotated  cone-shaped 
kernel  /).  The  SAD  may  be  written  as 

J-i 

Tk(i,  0;  zk)  =  2  CTj|vj  •  zk)i|2  (53) 

j=o 

where  Zk,i  is  the  (possibly  zero  extended)  spatial  sequence  surrounding  sensor  i  for  the  kth 
snapshot.  Due  to  the  rapid  drop  in  the  eigenvalues  of  commonly  used  kernels,  Tk(i,0;zk)  can 
be  well  approximated  as  T/t(i,0;Zk)  by  using  few  terms  J'  ( J '  <  J)  in  equation  (53).  This 
approximation  is  discussed  fully  in  [21]. 

Equation  (53)  can  be  interpreted  in  the  context  of  classical  beamforming  using  subarrays. 
A  subarray  beamformer  constructs  multiple  sub-beams  using  classical  beamforming  applied  to 
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subsets  of  sensor  elements  in  the  array.  Typically,  the  sensors  in  a  subarray  are  adjacent.  The 
sensors  corresponding  to  adjacent  subarrays  may  or  may  not  overlap  within  the  full  array.  A 
given  subarray  beamformer  estimates  the  spatial  spectrum  local  to  the  position  of  the  subarray 
within  the  full  array. 

At  each  sensor  index  i,  equation  (53)  is  a  weighted  sum  of  multiple  beamformers  evaluated 
using  the  same  spatial  signal  Zk,i  surrounding  sensor  i  (which  may  include  the  case  where  the 
full  spatial  signal  Zk  is  used).  This  weighted  sum  at  sensor  i  is  a  representation  of  the  spatial 
spectrum  at  sensor  i  (i.e.  the  SAD).  It  has  similar  form  to  a  classical  subarray  beamformer, 
although  differs  in  two  distinct  ways.  Up  to  J  beamformers  are  determined  at  each  sensor  (not 
one)  and  the  structure  of  the  beamformer  weights  and  the  combining  ratios  are  prescribed  by 
the  SAD  kernel  </>'.  It  is  important  to  note  that  the  J '  subarray  beamformers  can  be  applied  to 
Jl  consecutive  snapshots,  providing  that  the  source  maintains  its  signal  value  over  J'  samples. 
A  source  with  constant  modulus  property  will  naturally  satisfy  this  condition.  This  flexibility 
does  not  exist  in  a  typical  application  of  time-frequency  analysis.  A  deviation  from  the  constant- 
modulus  property  has  the  effect  of  changing  the  eigen-decomposition  by  scaling  the  eigenvalues 
in  equation  (53).  We  should  note  that  a  change  in  the  signal  phase  does  not  alter  the  magnitude 
squares  of  equation  (53).  The  ability  to  apply  each  term  of  the  approximation  in  equation  (53)  to 
a  different  snapshot,  instead  of  the  same  data  vector,  reduces  computational  time  and  hardware 
requirements. 


VI.  Simulations 

In  this  section,  we  demonstrate  the  proposed  approach  for  near-held  scatter  characterization. 
Consider  a  64  sensor  linear  equi-spaced  array  operating  at  a  carrier  frequency  of  6.14MHz  (A  = 
46.8m)  with  15m  spacing  between  sensors.  The  local  scatterer  distribution  comprises  a  point 
scatterer  in  the  near-held  at  range  1200m  and  angle  -30°,  from  the  array  center  (the  array  has 
total  length  945m).  Assume  the  test  source  angle  is  20°  with  respect  to  boresight  and  that 
the  test  source  is  temporally  stationary  and  coherent  with  and  26dB  stronger  power  than  the 
scattered  source.  In  this  example  we  have  used  the  Wigner  distribution  so  that  the  kernel  in  (12) 
is  5(m).  Of  course,  other  members  of  Cohen’s  class  may  also  be  used. 

Figure  4  shows  the  SAD  for  the  received  data.  The  plot  shows  the  received  power  (plot 
intensity)  as  a  function  of  angle  (spatial  frequency)  on  the  vertical  axis  and  sensor  position 
within  the  array  on  the  horizontal  axis.  The  spatial  frequency  axis  extends  beyond  the  interval 
±1  since  the  sensors  are  spaced  more  closely  than  The  SAD  is  dominated  by  the  substantially 
more  powerful  far-held  test  source  and  there  is  no  clear  indication  of  any  additional  scattering. 
The  far-held  source  has  the  same  angle  for  every  sensor,  and  therefore,  depicts  a  horizontal 
signature  in  the  s-a  domain.  In  hgure  5,  we  have  applied  the  orthogonal  projection  operator  to 
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the  received  sensor  data  and  computed  the  SAD  for  Pz^.  as  per  equations  (15)  and  (16).  The 
SAD  now  clearly  shows  the  presence  of  near-field  local  scatter.  In  this  case,  the  location  of  the 
near  field  source  is  far  enough  from  the  array  for  the  quadratic  phase  approximation  used  by 
Breed  and  Posch  to  be  valid  (i.e.  the  sensor  v.  angle  relationship  is  approximately  linear). 

The  beampatterns,  computed  using  a  70dB  Taylor  window,  for  the  cases  of  and  Pz^  are 
shown  in  figure  6.  Prom  the  plots,  the  presence  of  near-field  scatter  cannot  be  confirmed  as 
compared  with  alternative  explanations  for  the  distorted  beampatterns,  such  as  poor  array  cali¬ 
bration. 

Now  we  wish  to  demonstrate  the  use  of  the  alias-free  summed  beamformer  SAD  described 
in  equation  (53).  Again  the  simulations  assume  a  64  element  linear  equi-spaced  array,  with 
operating  frequency  and  sensor  spacing  as  before. 

In  the  first  simulation,  a  far-field  source  was  placed  at  angle  30°  (all  angles  are  with  respect 
to  array  boresight)  and  a  near-field  source  placed  at  a  range  of  100m  and  an  angle  of  -40°.  The 
Choi- Williams  SAD  (a  =100)  implemented  using  the  alias-free  rotated  kernel  form  is  shown  in 
figure  7.  Signals  apparently  arriving  from  spatial  frequencies  greater  than  1  or  less  than  -1  do 
not  correspond  to  real  signals.  The  crosses  marked  on  the  figure  indicate  the  true  angles  of  each 
source  at  every  sensor. 

The  weighted  summed  beamformer  implementation  of  the  Choi-Williams  SAD  described  in 
equation  (53)  has  been  computed  for  the  same  scenario  as  in  figure  7.  The  result  is  shown  in 
figure  8  for  the  case  where  J'  =25  and  J  =65  (i.e.  25  beamformers  have  been  retained  out  of  a 
possible  total  of  65).  Note  that  visually,  there  is  no  apparent  loss  of  fidelity  using  the  approximate 
SAD.  This  implementation  has  substantially  lower  computational  complexity. 

The  simulation  results  shown  in  figures  9  and  10  are  for  the  case  of  two  near-field  sources.  In 
figure  9,  one  source  is  placed  at  a  range  of  300m  and  an  angle  of  30°,  whereas  the  second  is  at 
a  range  of  200m  and  an  angle  of  -40°.  The  SAD  is  a  Choi-Williams  kernel  implemented  using 
weighted  summed  beamforming  with  J'=  25  and  J  =65  (i.e.  25  of  65  retained).  Note  the  presence 
of  cross  terms  between  the  two  source  sensor-angle  signatures.  The  Choi-Williams  kernel  is  well 
known  to  perform  poorly  in  terms  of  cross  term  suppression  for  the  case  when  the  underlying 
signals  are  other  than  close  to  parallel  to  either  the  sensor  (time)  axis  or  the  angle  (frequency) 
axis  in  the  SAD  (TFD). 

The  results  presented  in  figure  10  again  show  the  result  of  a  SAD  using  a  Choi-Williams  kernel 
implemented  with  the  25  term  weighted  summed  beamformer.  There  are  two  near-field  sources 
and  the  source  positions  have  been  selected  to  be  along  the  axis  of  the  array.  The  first  source 
is  at  300m  and  90°  while  the  second  is  at  600m  and  -90°.  The  first  source  is  embedded  within 
the  array  and  co-located  exactly  at  a  sensor  position.  The  second  source  is  beyond  the  lefthand 


15 


extent  of  the  array.  The  first  source  has  been  localized  at  exactly  the  correct  sensor  within 
the  array  (at  300m  from  the  array  center).  Near-field  spreading  loss  (1/r)  between  sensors  has 
caused  the  source  intensity  to  drop  dramatically  at  adjacent  sensors.  The  corresponding  dual 
time-frequency  domain  signature  is  that  of  an  impulse.  The  second  source  is  again  correctly 
represented  as  being  at  -90°  to  the  left  end  of  the  array.  Notice  that  the  1/r  spreading  loss  for 
the  near-field  source  means  that  the  SAD  has  stronger  intensity  at  the  left  end  of  the  array.  This 
is  the  degenerate  near-field  situation,  where  even  though  the  source  is  in  the  near-field  of  the 
array,  the  angle  to  each  sensor  in  the  array  is  identical. 

The  simulations  shown  in  figures  11  through  16  contrast  the  cases  of  point  sources  in  the  far 
and  near-field  and  a  spatially  spread  source  in  the  far-field.  Traditional  beampatterns  (computed 
using  a  70dB  Taylor  weight)  show  that  it  may  be  difficult  to  discriminate  between  spatially 
spread  far-field  sources  and  a  near-field  point  source.  Using  the  SAD,  one  is  clearly  able  to 
visually  separate  these  cases  (contrast  figures  14  and  16).  Note  that  the  spatially  spread  source 
was  generated  using  20  realizations  of  a  low  pass  noise  process  modulating  the  spatial  steering 
vector  corresponding  to  the  mean  angle  for  the  spread  source. 

VII.  Conclusions 

We  have  proposed  a  new  method  for  characterizing  near-field  scatter  local  to  a  receiving  array. 
As  part  of  the  characterization,  we  have  exploited  and  generalised  the  spatial  Wigner  distribution 
proposed  by  Breed  and  Posch,  although  we  have  renamed  it  the  sensor  angle  distribution  to  avoid 
confusion  with  a  similarly  named  but  differently  defined  spatial  time-frequency  distribution.  An 
orthogonal  projection  operator  derived  from  the  steering  vector  for  the  far-field  test  source  is 
used  to  exclude  the  direct  propagation  path  from  the  test  source  in  the  characterization. 

This  paper  has  addressed  several  issues  concerning  the  use  of  the  Sensor-Angle  Distribution 
(SAD)  as  a  tool  for  characterizing  near-field  scatter  environments. 

An  expression  has  been  derived  describing  the  angle  to  a  source  or  scatter  site  from  each 
sensor  in  an  array.  Formulas  are  given  for  least  squares,  and,  under  certain  conditions,  maximum 
likelihood  estimates  of  source  location  based  on  source  angle  measurements  such  as  those  obtained 
from  a  SAD. 

The  use  of  a  rotated  smoothing  kernel  has  proven  effective  at  overcoming  aliasing  in  the  SAD 
caused  by  the  spatial  frequencies  obtained  for  realistically  positioned  sources. 

An  approximate  method  for  implementing  the  SAD  has  been  found  to  have  a  direct  interpre¬ 
tation  in  terms  of  a  multiple  subarray  beamformer.  However,  the  weighting  functions  applied  to 
the  beamformer,  and  the  fact  that  more  than  one  beamformer  is  combined  to  construct  the  final 
spatial  spectrum  at  each  sensor  is  an  extension  of  present  subarray  concepts. 

The  SAD  presented  in  this  paper  assumed  the  use  of  a  data-independent  sensor-angle  kernel, 
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as  reflested  in  equation  (12).  It  is  important  to  note  that  data-dependent  kernels  can  also  be 
employed  in  the  SAD.  Positive  time- frequency  distributions  [22],  optimum  time- frequency  kernels 
[23],  and  the  time- frequency  reassignment  method  [24]  can  all  be  cast  in  terms  of  the  two  new 
variables,  sensor  and  angle,  leading  to  potential  improvements  over  equation  (12)  in  resolution 
and  cross-term  suppression. 
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Appendix 

I.  Cramer-Rao  Bound  for  Location  from  Sensor  Measurements 

Consider  the  received  signal  model 

zk  =  a  •  sk  +  wk  (54) 

where  k  is  the  array  snapshot  index.  For  the  case  of  a  single  near-field  source  at  (rs,0s)  then  sk 
is  a  scalar  s k  and  a  is  given  by  equation  (8).  Let  the  total  number  of  snapshots  be  N  with  the 
index  k  E  {1, . . .  ,  N}.  The  noise  is  additive  Gaussian  noise  where  wk  =  [wi, . . .  ,wm]J  which 
for  all  received  snapshots  can  be  written  W  =  [wi,...  ,wn],  The  noise  is  characterized  by 
E[wijkWjc/jk/]  =  <7?k^i-i'5k-k'  =  of,  or,  more  completely  ^  ~  Af(  0,  of).  In  this  case  the  variance 
of  ith  sensor  may  be  different  for  each  sensor  i,  but  remains  the  same  over  all  snapshots  for  a  given 
sensor.  Let  Z  =  [zi, . . .  ,  zn]  be  the  data  matrix  of  all  received  data.  Let  s  =  [si, . . .  ,  sn]  be  the 
deterministic  source  signal  characterized  only  by  the  vector  s.  Equation  (54)  can  be  re-written 
as 


Z  ~  a(rs,  0S)  •  s  +  W  (55) 

which  has  unknown  parameters  (rs,05)  and  nuisance  parameter  s  (which  increases  in  dimension 
as  N  increases). 

The  likelihood  function  for  this  model  is 

N 

p(Z ;  rs,  0S,  s)  =  p(zk  5  rsj  sk)  (56) 

k=l 

with 

^  1  /  i  \ 

p(zk  ;rs,es, sk)  =  TT  - exp  -—2  [zi>k  -  ai(rs,0s)  •  sk]2  )  (57) 

i=i  pita?  V  2(7.  J 
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hence 


N  M  1  /  1  \ 
p(  Z;rs,0s,s)  =  Yl 


(58) 


with  the  log-likelihood  given  by 


M 


M  r 


2^.2  [Zi»k  ai(rs?^s)  *  sk] 


(59) 


lnp(Z ;  rs,  6S ,  s)  =  -N  ^  ln[-v/27rcri]  -  N  ^ 

i=l  i=l  L  i 

The  FIM  is  constructed  from  the  partial  derivatives  of  ai(rs,6s )  •  sfc  with  respect  to  rs,  9S  and 
sk.  Let 


or 


Ui,k  —  si  @s)  '  $k 


Sk  —j22Lr  * 

Vi,k  = - e  1  x  r*.« 

rs,i 


Inserting  the  r;)5  gives 

Vi  ,h  = 


Sk 


x2  —  2  rsa:  sin  9S  +  r2 
where  x  at  each  sensor  position  is 

M  —  1 


.  e-iy\/x!-2rJisin«s  +rf 


X  =  d- 


(®  1) 


for  *  G  ,A/} 


The  partial  derivatives  are  (in  terms  of  x  to  improve  clarity) 
dyi,k  _  sk  ■  e-jx\/l2-2r^sin^+rs2 


drs 


x2  —  2  rsx  sin  6S  +  r2 


[-xsin0s  +  r5] 


-1 


.2  7T 


_ _ _  —  J  - 

\Jx2  —  2  rsx  sin  9S  +  r2  A 


(60) 


(61) 


(62) 


(63) 


(64) 


dVi,k  _  Sk  ■  e-JxV^r^sin^+rf 
<9#s  ar2  —  2  rsx  sin  9S  +  r2  C°S  s 


&yi,k 


^Jx2  —  2ria;sin0,s  +  r2 


-  J  -y  y/  x2-2  /yx  sin  $s4-r2 


+  3- 


,2n 


(65) 


dsk  ijx2  —  2rsx  sin  9S  +  r2 
For  the  case  of  A:  =  1  (i.e.  the  single  snapshot  case),  the  individual  elements  of  the  FIM  are 


(66) 


I(rsA,sk)  = 


Em  i  dyjjc  dyjjc 
2=1  ~af  drs  8rs 

EM  l  dy^k  dyj,k 
i=l  aj  89 s  drs 

EM  1  dyitk  dyj,k 
i=l  of  8sk  drs 


Em  1  8yi  k  8yitk  v^M  1  8yi,k  dyitk 

2—1  - _  ^ - “ 


t?  8rs 


EM  X  dyiik  dyitk 

2=1  ^  d6s  80s 

M  1  dyitk  8yiik 


Em  i 
2=1  ^ 


of  8s  k  89 s 


2-/2=l  ^f  dr s  8sk 

Em  1  dyj^k  8yj  k 

2=1  89s  8sk 

EM  l  8yiyk  dyiik 
i=l  of  8sk  8sk 


(67) 
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The  CRLB  for  each  of  the  unknown  parameters  is  determined  from  the  inverse  of  equation  (67) 
evaluated  using  equations  (64)  through  (66)  at  the  true  values  of  the  unknown  parameter. 


var(rs)  >  p(r„0s,sk)]lti 

(68) 

var(0s)  >  [I(rs,0s,Sk)]2)2 

(69) 

var(sfc)  >  [I(rs,0s,sk)]3  3 

(70) 

Appendix 

I.  Cramer-Rao  Bound  for  Location  from  Sensor  Phase  Measurements 

Consider  the  case  where  only  the  phase  of  the  received  signal  is  recorded  at  each  sensor.  In 
this  case  equation  (61)  becomes 


Ui}k 


=  e  3  a  Ts 


(71) 


where,  as  before,  substituting  in  r^s  gives 

yik  =  e-if^2-2rsMn9s+r2  (72) 

The  partial  derivatives  of  the  two  remaining  unknown  parameters  (rs,0s)  are  (in  terms  of  sensor 
location  x) 

dyi,k  _  _  .27r 
drs  ~  3  A 


-xsinQs  +  rs 


W x2  —  2rsx  sin05  +  rs2 


-  j  If  Vx2  ~ 2  rs  x  sin  Os  +rf 


(73) 


dm*  _  .27T 
3  A 


rs  x  cos  0S 


_Vz2  —  2r5  £sin0s  +  rs2\ 

The  elements  of  the  FIM  are  then 

„  ,  _V  1  ‘,ykk  ',ykk 

I— 1  1 


j~  Vx2~2rs  x sin  03+r% 


e  J  * 


(74) 


(75) 


Tfr  fl  1  V  1  ^ 

I(rs,0.)i,2-2^Z2  ST"3T 


.  1  a?  drs  d0s 

1=1  1 


(76) 


M 


I(rs,0s)2,i  =  ^ 


l  3yi,kdyi,k 


i=i  » 


a?  80 *  <9rs 


(77) 
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(78) 


T(r  0)  =  V  1  dyj.k  dyi>k 

(  S’0s)2,2  Z^n-2 


.  ,  a?  des  dos 

1=1  1 


The  Cramer-Rao  lower  bounds  for  estimates  of  rs  and  0S  based  on  phase  only  measurements  at 
each  sensor  in  the  array,  are,  respectively 

_ I(r„fl,)2,a _ 


and 


var(fs)  > 


var (6S)  > 


[I(rs,  0s)i,i  •  I(rs, @s)2,2  I(r„ 08)i,2  ■  I(i*S) $s)2,i] 

_ I(rs,  fl8)i,i _ 

[I(rs,  0s)i,i  •  I(fs)  @5)2,2  -  I(r8, 0s)i,2  •  I(rs,  #s)2,i] 


(79) 


(80) 
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Fig.  1.  Sensor  array  geometry  and  notations  for  a  linear  array  and  a  near-field  source 
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Fig.  2.  CRLB  (in  m2)  of  the  range  to  a  source,  for  various  source  locations,  for  the  case  of  a  64  sensor 
linear  array,  spaced  15m  between  sensors.  The  sensor-angle  is  assumed  to  have  been  estimated  at 
each  sensor  with  a  variance  of  1  deg2. 
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Fig.  3.  CRLB  (in  deg2)  of  the  angle  to  a  source,  for  various  source  locations,  for  the  case  of  a  64  sensor 
linear  array,  spaced  15m  between  sensors.  The  sensor-angle  is  assumed  to  have  been  estimated  at 
each  sensor  with  a  variance  of  1  deg2. 


sensor  position  wit  array  center  (m) 

Fig.  4.  SAD  for  the  received  sensor  data  Zk.  The  far-field  test  source  (at  angle  20°)  dominates  the  SAD 
characterization. 


sensor  position  wit  array  center  (m) 

Fig.  5.  SAD  for  the  received  data  following  application  of  the  orthogonal  projection  operator,  Pzk. 
With  the  direct  propagation  path  from  the  far-field  test  signal  removed  by  the  orthogonal  projection 
operator,  the  local  scatterer  spectrum  is  revealed.  The  local  point  scatterer  is  at  a  range  of  1200m 
and  angle  of  -30°. 
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Beampattern 


Fig.  6.  Beampattern  for  zk  (-  )  and  Pzk  ( — ).  Without  the  sensor  angle  characterization  it  is  not 
possible  to  identify  perturbations  from  the  ideal  test  source  beampattern  as  being  due  to  near-field 
scatter.  For  example,  poor  array  calibration  may  be  indistinguishable  from  local  scatter. 


One  FF  and  One  NF  source  (CW  sigma=tOO) 


Fig.  7.  Choi- Williams  SAD  for  a  far-held  source  at  angle  30°  and  a  near-field  source  at  range  100m  and 
angle  -40°.  True  angles  for  each  source  at  every  sensor  are  marked  by  crosses. 
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One  FF  and  One  NF  source  (CW  sigma- 100) 


Fig.  8.  Approximate  weighted  summed  beamformer  SAD  for  a  far-field  source  at  angle  30°  and  a  near¬ 
field  source  at  range  100m  and  angle  -40°.  In  this  case  25  of  a  possible  65  beamformers  are  summed 
together  (see  equation  (53)). 


Two  NF  sources  (CW  sigma=100) 


sensor  position  wit  array  center  (m) 

Fig.  9.  Approximate  weighted  summed  beamformer  results  for  two  near-field  sources  at  300m  and  30° 
and  200m  and  -40°.  25  of  a  possible  65  beamformers  are  summed  together. 
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Two  NF  sources  (CW  sigma=100) 


sensor  position  wrt  array  center  (m) 


Fig.  10.  Approximate  weighted  summed  beamformer  results  for  two  near-field  sources  where  the  sources 
are  aligned  along  the  axis  of  the  array  at  600m  and  -90°  and  300m  and  90°. 
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Fig.  11.  Beampattern  for  a  far-field  point  source  at  angle  20°. 
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One  far -field  point  source  (CW  sigma=100) 


Fig.  12.  Approximate  weighted  summed  beamformer  SAD  (J'=25  and  J=  65)  result  for  a  far-field  point 
source  at  angle  20°. 


One  near-field  point  source  (Beampattem) 


Fig.  13.  Beampattem  for  a  near-field  point  source  at  range  800m  and  angle  20°. 
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sensor  positron  wit  srrsy  center  (m) 


Fig.  14.  Approximate  weighted  summed  beamformer  SAD  (25  of  65)  result  for  a  near-field  point  source 
at  range  800m  and  angle  20°. 


Fig.  15.  Beampattern  for  a  far-field  spatially  spread  source  with  mean  angle  20°. 
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One  far-field  spread  source  (CW  sJgma=100> 
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Fig.  16.  Approximate  weighted  summed  beamformer  SAD  (25  of  65)  result  for  a  far-field  spatially  spread 
source  with  mean  angle  20°. 
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Improved  Blind  Separations  of 
Nonstationary  Sources  Based  on 
Spatial  Time-Frequency  Distributions 


Yimin  Zhang  and  Moeness  G.  Amin 


Abstract 

Blind  source  separation  based  on  spatial  time-frequency  distributions  (STFDs)  provides  improved 
performance  over  blind  source  separation  methods  based  on  second-order  statistics,  when  dealing  with 
nonstationary  signals  that  are  localizable  in  the  time-frequency  (t-f)  domain.  In  the  existing  STFD  meth¬ 
ods,  the  covariance  matrix  is  first  used  to  whiten  the  data  vector,  then  the  mixing  matrix  and  subsequently 
the  source  waveforms  are  estimated  using  STFD  matrices  constructed  from  the  source  t-f  autoterm  points. 
This  letter  improves  the  STFD-based  source  separation  method  by  performing  both  whitening  and  esti¬ 
mation  steps  using  the  source  t-f  signatures.  This  modification  provides  robust  performance  to  noise,  and 
allows  reduction  of  the  number  of  sources  considered  for  separation. 
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I.  Introduction 


Several  narrowband  blind  source  separation  (BSS)  methods  have  been  proposed  in  the 
literature  [1],  [2],  [3],  [4].  Generally,  blind  source  separations  for  independent  sources  are 
performed  based  on  the  employment  of  at  least  two  different  sets  of  matrices  that  span  the 
same  signal  subspace.  One  matrix  is  used  for  whitening  purpose,  while  others  are  jointly 
used  to  estimate  the  spatial  signatures  and  source  waveforms  impinging  on  a  multi-antenna 
receiver.  Covariance  matrices  with  different  time  lags,  or  cumulants  matrices  with  different 
orders  are  typically  used  for  the  above  purpose. 

Nonstationary  signals  are  frequently  encountered  in  radar,  sonar,  and  acoustic  applica¬ 
tions  [5],  [6].  For  nonstationary  signals,  time-frequency  distributions  (TFDs)  have  been 
recently  employed  for  array  processing  and  found  successful  in  blind  separations  of  nonsta¬ 
tionary  signals  [7],  [8],  [9],  [10],  [11],  [12].  These  methods  are  particularly  effective  when 
the  signals  are  highly  localized  in  the  time-frequency  (t-f)  domain.  In  these  methods,  while 
the  spatial  time-frequency  distribution  (STFD)  matrices  are  used  for  source  diagonaliza- 
tion  and  anti-diagonalization,  the  whitening  matrix  remains  the  signal  covariance  matrix. 
The  STFD  matrices  are  constructed  from  the  auto-TFDs  and  cross-TFDs  of  the  sensor 
data  and  evaluated  at  different  points  of  high  signal-to-noise  ratio  (SNR)  pertaining  to  the 
t-f  signatures  of  the  sources.  Although  this  method  improves  the  performance  over  that 
based  on  solely  the  covariance  matrices,  yet  it  does  not  fully  utilize  all  inherent  advantages 
of  STFD. 

This  letter  modifies  the  existing  STFD-based  BSS  methods  by  performing  whitening 
using  an  averaged  STFD  matrix  over  the  source  t-f  signatures  and,  as  such,  making  use  of 
the  t-f  localization  properties  of  the  sources  in  both  the  whitening  step  and  joint  estimation 
step  of  the  source  separation  procedure. 

Employing  the  STFD  for  whitening  leads  to  robustness  of  subspace  decompositions  to 
noise  and,  thereby,  enhances  the  unitary  mixture  representations  of  the  problem.  It  also 
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applies  equal  source  discrimination  and  allows  the  consideration  of  the  same  subset  of 
sources  in  both  the  whitening  and  joint-diagonalization  phases. 

II.  Signal  Model 

In  narrowband  array  processing,  when  n  signals  arrive  at  an  m-e lement  array,  the  linear 
data  model 

x(i)  =  y(t)  +  n(i)  =  Ad(t)  +  n (t)  (1) 

is  commonly  used,  where  A  is  the  mixing  matrix  of  dimension  mxn  and  is  assumed  to  be 
full  column  rank,  x(t)  =  [xi(t),  ...,xm(t)]T  is  the  sensor  array  output  vector,  and  d (t)  = 
[di(t), ...,  dn(t)]T  is  the  source  signal  vector,  where  the  superscript  T  denotes  the  transpose 
operator.  n(t)  is  an  additive  noise  vector  whose  elements  are  modeled  as  stationary, 
spatially  and  temporally  white,  zero-mean  complex  random  processes,  independent  of  the 
source  signals. 

The  covariance  matrix  is  given  by 

Rxx  =  i?[x(i)xH(t)]  =  Ryy  +  al  =  ARddA^  +  al,  (2) 

where  E(-)  is  the  statistical  expectation  operator,  the  superscript  H  denotes  conjugate 
transpose,  Rdd  —  £'[d(t)d//(t)]  is  the  signal  covariance  matrix,  a  is  the  noise  power  at 
each  sensor,  and  I  denotes  the  identity  matrix.  We  assume  that  Rxx  is  nonsingular,  and 
the  observation  period  consists  of  N  snapshots  with  N  >  m. 

III.  Blind  Source  Separation  based  on  Spatial  Time-Frequency  Signatures 

A.  Spatial  Time- Frequency  Distributions 

The  discrete  form  of  Cohen’s  class  of  STFD  of  the  data  snapshot  vector  x(f)  is  given 

by  [7] 

°o  CO 

D„(<,/)=  £  £  ^,T)x((  +  i  +  T)x"((  +  i-r)e-^'  (3) 

<=—00  T—  OO 
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where  <f>(l,r)  is  a  t-f  kernel.  Substituting  (1)  into  (3),  we  obtain 


Dxx(^)  /)  -  Dyy  (t,  /)  +  Dyn(t,  /)  +  Dny(t,  /)  +  DlUl(U  /)•  (4) 

Under  the  uncorrelated  signal  and  noise  assumption  and  the  zero-mean  noise  property, 
E  [Dyn(t,  /)]  =  E  [D„y(t,  /)]  =  0.  It  follows 

E  [Dxx(i,  /)]  =  Dyy(t,  f)  +  E  [Dnn(f ,  /)]  =  ADdd(t,  f)AH  +  E  [Dnn(t,  /)] .  (5) 

Similar  to  (2),  which  relates  the  signal  covariance  matrix  to  the  data  spatial  covariance 
matrix,  (5)  provides  the  basis  for  source  separation  by  relating  the  STFD  matrix  to  the 
source  TFD  matrix,  Ddd(t,  /). 

B.  Blind  Source  Separation  [7] 

In  blind  source  separation  techniques,  there  is  an  ambiguity  with  respect  to  the  order 
and  the  complex  amplitude  of  the  sources.  It  is  convenient  to  assume  that  each  source  has 
unit  norm,  that  is,  Rdd  =  I. 

The  first  step  of  blind  source  separation  based  on  TFDs  is  whitening  of  the  signal  x(t)  of 
the  observation.  This  is  achieved  by  estimating  the  noise  power  1  and  applying  a  whitening 
matrix  W  to  x(t),  i.e.,  an  n  x  m  matrix  satisfying: 

WRyyWH  =  W(RxX  -  crl) WH  =  WAAhWh  =  I.  (6) 

The  whitening  matrix  is  estimated  from  the  eigendecomposition  of  Rxx  [7].  The  accuracy 
of  the  whitening  matrix  estimate  depends  on  the  estimation  accuracy  of  the  eigenvectors 
and  eigenvalues  corresponding  to  the  signal  subspace.  The  whitened  process  z (t)  —  Wx(t) 
still  obeys  a  linear  model, 

z  (t)  =  Wx(t)  =  WAd(t)  +  Wn .(f)  =  Ud(t)  +  Wn  (t),  (7) 

lrThe  noise  power  can  be  estimated  only  when  m  >  n  [7].  If  m  =  n,  the  estimation  of  the  noise  power  becomes 
unavailable  and  a  —  0  will  be  assumed. 
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where  U  A  WA  is  an  n  x  n  unitary  matrix. 

The  next  step  is  to  estimate  the  unitary  matrix  U.  The  whitened  STFD  matrices  in 
the  noise-free  case  can  be  written  as 


D  „(t,  f)  =  WDxx(t,  f)WH  =  UD  Ad(t,  f)XJH.  (8) 

In  the  autoterm  regions 2,  Ddd(i,  /)  is  diagonal,  and  an  estimate  U  of  the  unitary  matrix  U 
may  be  obtained  as  a  joint  diagonalizer  of  the  set  of  whitened  STFD  matrices  evaluated  at 
K  autoterm  t-f  points,  {D22(tj,  fi)\i  =  1, ...,  K}.  The  source  signals  and  the  mixing  matrix 
can  be,  respectively,  estimated  as  d  (t)  =  UWx(t)  and  A  =  W#U,  where  superscript  * 
denotes  pseudo-inverse. 

IV.  Properties  of  STFD  Matrices 

In  reference  [13],  the  subspace  analyses  of  STFD  matrices  are  presented  for  signals  with 
clear  t-f  signatures,  such  as  frequency  modulated  (FM)  signals.  It  was  shown  that  the 
offerings  of  using  a  STFD  matrix  instead  of  the  covariance  matrix  are  two-fold.  First,  the 
selection  of  autoterm  t-f  points,  e.g.,  points  on  the  source  instantaneous  frequency  where 
the  power  is  concentrated,  enhances  the  SNR.  Second,  the  difference  in  the  t-f  localization 
properties  of  the  source  signals  permits  source  discrimination  and  allows  the  selection  of 
fewer  sources  for  matrix  construction.  In  the  presence  of  noise,  the  consideration  of  a 
subset  of  signal  arrivals  reduces  perturbation  in  matrix  eigendecomposition,  and  becomes 
essential  when  there  is  insufficient  number  of  sensors. 

The  prime  motivation  of  this  letter  is  to  make  use  of  the  above  properties  in  the  whiten¬ 
ing  phase  of  TFD-based  BSS  methods.  In  particular,  the  implicit  reduction  in  the  noise 
as  a  result  of  autoterm  selection  and  fewer  source  considerations  lead  to  improved  signal 
subspace  and  source  number  estimation.  Both,  according  to  (11),  are  key  to  the  whitening 
matrix  construction. 

2The  selection  of  autoterm  regions  has  been  discussed  in  [11],  [12],  [14]. 
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V.  Modified  Method  for  Source  Separation 

In  the  method  proposed  in  [7]  and  summarized  in  Section  III,  although  STFD  matrices 
are  used  to  estimate  the  unitary  matrix  U,  the  covariance  matrix  is  still  used  in  the 
whitening  process.  Using  the  STFD  matrix  in  place  of  the  covariance  matrix  Rxx  to 
perform  whitening  is  an  attractive  alternative. 

Denote  Dxx(ti,  /i), ...,  Dxx(i/f,  fx)  as  the  STFD  matrices  constructed  from  K  autoterm 
points  being  defined  over  a  t-f  region  Q  and  belonging  to  fewer  n0  <  n  signals.  Also, 
denote,  respectively,  d°(t)  and  d  (t)  as  the  nQ  and  n  —  n0  sources  being  present  and  absent 
in  the  t-f  region  Q.  The  n  —  n0  sources  could  be  undesired  emitters  or  sources  to  be 
separated  in  the  next  round  of  processing.  The  value  of  n0  is  generally  unknown  and  can 
be  determined  from  the  eigenstructure  of  the  STFD  matrix.  Using  the  above  notations, 
we  have 

x(t)  =  A°d°(t)  4-  Ad  (t)  +  n(t),  (9) 

where  A°  and  A  are  the  mxn0  and  m  x  (n  —  n0)  mixing  matrices  corresponding  to  d°(f) 
and  d  (t)  ,  respectively. 

Let  Dxx  be  the  average  STFD  matrix  of  a  set  of  STFD  matrices  defined  over  the 
same  region  D,  but  may  incorporate  a  different  t-f  kernel.  The  incorporation  of  multiple 
(t,  f )  points  through  the  averaging  process  reduces  the  noise  effect  on  the  signal  subspace 
estimation,  as  discussed  in  [13].  Denote  atf  as  the  estimation  of  the  noise-level  eigenvalue 
of  the  STFD  matrix  Dxx.  Then, 

WDyyWfl  =  W(DXX  -  atfI)WH  =  WA°Djd  (WAf  =  I.  (10) 

In  (10),  due  to  the  ambiguity  of  signal  complex  amplitude  in  BSS,  we  have  assumed  that 
the  averaged  source  TFD  matrix  D£d  corresponding  to  d °(t)  is  I  of  n0  x  n0  for  convenience 
and  without  loss  generality.  Therefore,  the  whitening  matrix  W  is  obtained  as 

W  =  [(Af  -  ^)-,/2 hi', ....  (A l  -  ” ,  (11) 
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where  A^, are  the  n0  largest  eigenvalues  of  Dxx  and  h4/, hlJo  are  the  corresponding 
eigenvectors  of  Dxx.  Note  that  Ddd  and  Dyy  are  of  reduced  rank  n0  instead  of  rank  n,  due 
to  the  source  discrimination  performed  through  the  selection  of  the  t-f  points  or  specific 
t-f  region.  Therefore,  WA°  =  U  is  a  unitary  matrix,  whose  dimension  is  n0  x  na  rather 
than  nx  n.  The  whitened  process  z(t)  becomes 

z(t)  =  Wx(t)  =  WA°d°(t)  +  WAd  (t)  +  Wn  (t)  =  Ud°(t)  +  WAd(f)  +  Wn{t),  (12) 

♦ 

In  the  t-f  region  Cl,  the  TFD  of  d (t)  is  zero  and,  therefore,  the  averaged  STFD  matrix  of 
the  noise-free  components  becomes  an  identity  matrix,  i.e., 

Dzz  =  WDXXW"  =  UD2dUH  =  I.  (13) 

Eqn.  (13)  implies  that  the  autoterm  and  crossterm  TFDs  averaged  over  the  t-f  region 
Cl  become  unity  and  zero,  respectively,  upon  whitening  with  matrix  W.  U  as  well  as 
the  mixing  matrix  and  source  waveforms  are  estimated  following  the  same  procedure  of 
Section  III.  It  is  noted  that,  when  na  =  1,  source  separation  is  no  longer  necessary. 

Selection  of  the  same  number  of  sources,  n0,  should  be  done  at  both  whitening  and  joint 
diagonalization  stages,  otherwise  mismatching  of  the  corresponding  sources  will  result. 
While  our  proposed  modified  blind  source  separation  method  provides  the  mechanism  to 
satisfy  this  condition,  the  covariance-based  whitening  approach  does  not  lend  itself  to 
avoid  the  mismatching. 

VI.  Simulation  Results 

We  consider  a  two-sensor  array  with  a  half-wavelength  spacing.  For  simplicity,  we 
consider  chirp  signals  as  the  sources  to  be  separated.  In  the  first  example,  two  chirp 
signals  are  received  in  the  presence  of  white  Gaussian  noise.  The  sources  arrive  from 
different  directions  9\  =  30°  and  92  =  35°.  The  normalized  start  and  end  frequencies  of 
the  first  chirp  signal  are  0  and  0.5,  respectively,  whereas  those  for  the  second  signal  are 
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0.1  and  0.4.  The  number  of  snapshots  used  is  1024. 

The  performance  is  evaluated  by  using  the  mean  rejection  level  (MRL),  defined  as  [7] 

MRL  =  ££|(A#A)P9|2  (14) 

where  A  is  the  estimate  of  A.  A  smaller  value  of  the  MRL  implies  better  source  separation 
results.  An  MRL  lower  than  —10  dB  is  considered  satisfactory  [7]. 

Fig.  1  shows  the  MRL  versus  the  input  SNR  of  the  two  sources  (we  assume  that  all  signals 
have  the  same  power).  The  dashed  line  corresponds  to  the  existing  method  where  the 
covariance  matrix  Rxx  is  used  for  whitening,  and  the  solid  line  corresponds  to  the  modified 
method  where  the  averaged  STFD  matrix  Dxx  is  used  instead.  In  the  latter,  the  average 
of  spatial  pseudo-Wigner-Ville  distributions  (SPWVDs)  of  window  size  257  is  applied  to 
estimate  the  whitening  matrix.  For  the  estimation  of  the  unitary  matrix  for  both  methods, 
the  spatial  Wigner-Ville  distribution  (SWVD)  3  matrices  using  the  entire  data  record  are 
computed.  The  number  of  points  used  to  perform  the  joint  diagonalization  is  K  =  32  for 
each  signal  and  the  points  are  selected  on  the  signal  instantaneous  frequencies.  The  curves 
are  calculated  by  averaging  100  independent  trials  with  different  noise  sequences.  Fig.  1 
clearly  shows  the  improvement  when  STFDs  are  used  in  both  phases  of  source  separations, 
specifically  for  low  SNRs.  To  satisfy  the  —10  dB  MRL,  the  required  input  SNR  is  about 
17.8  dB  for  the  existing  method  and  is  about  5.6  dB  for  the  modified  method. 

To  examine  the  advantages  of  source  discrimination  in  the  t-f  domain,  we  add  another 
chirp  signal  to  the  above  scenario.  The  third  chirp  arrives  from  direction  63  =  40°,  and 
its  normalized  start  and  end  frequencies  are  0.15  and  0.55,  respectively.  Other  parameter 
settings  are  the  same  as  those  used  in  Fig.  1.  Only  the  t-f  points  belonging  to  the  first 
two  signals  are  chosen  to  construct  the  STFD  matrices  for  both  methods  (i.e.,  n0  =  2). 
It  is  evident  from  Fig.  2  that,  by  discriminating  against,  or  filtering  out,  the  third  signal, 

3Tbe  method  proposed  here  is  not  limited  to  use  specific  TFDs  and  the  SPWVD  and  SWVD  are  chosen  for 
simplicity.  Other  TFDs  can  also  be  used. 
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the  modified  method  provides  very  close  performance  to  the  two-signal  case,  whereas  the 
existing  covariance  matrix  based  whitening  method,  due  to  the  inability  to  discriminate 
among  sources  and  provide  sufficient  degrees-of-freedom  for  subspace  estimation,  cannot 
provide  satisfactory  MRL,  even  when  the  input  SNR  is  high. 

VII.  Conclusion 

A  modified  blind  source  separation  method  for  nonstationary  signals  is  introduced. 
In  this  method,  the  TFD  signal  localization  properties  are  fully  utilized  for  improved 
separation  performance.  The  proposed  modification  provides  noise  robustness  and  can  be 
used  to  reduce  the  number  of  sources  considered  for  separation. 
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and  Cyclostationary  Signals 
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Abstract 

The  cyclic  autocorrelation  function,  commonly  used  for  cyclostationary  signals,  and  the  ambiguity 
function,  typically  employed  for  analysis  and  recovery  of  nonstationary  signals,  such  as  FM,  have  the 
same  formulation.  However,  nonstationary  and  cyclostationary  signals  have  distinct  localization  proper¬ 
ties  in  the  time-lag  frequency-lag  domain.  Therefore,  nonstationary  and  cyclostationary  signals  can  be 
represented  and  processed  within  the  same  framework  for  many  applications,  with  the  distinct  signatures 
allowing  effective  source  discriminations.  An  example  in  array  processing  is  given  where  nonstationary 
and  cyclostationary  signals  are  separated  following  simple  spatial  signature  estimation  exploiting  the 
aforementioned  properties. 
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I.  Introduction 


Most  digital  communication  signals  are  stationary  when  they  are  sampled  at  the  symbol- 
rate.  On  the  other  hand,  when  they  are  fractionally  sampled,  they  demonstrate  the  so- 
called  cyclostationary  property,  that  is,  the  self-correlation  function  of  such  a  signal  is 
periodic  with  non-zero  cyclic  frequencies  [1].  Exploitation  of  the  cyclostationarity  allows 
discrimination  of  source  signals  if  their  cyclic  frequencies  are  distinct.  This  property 
has  been  widely  discussed  for  various  purposes  in  array  processing,  for  example,  source 
separation  and  direction  finding  [2],  [3],  [4],  [5]. 

Another  function  with  similar  definition  is  the  ambiguity  function  [6].  An  ambiguity 
function  is  usually  used  for  nonstationary  signal  analysis  and  waveform  design.  The  con¬ 
sideration  of  the  ambiguity  function  for  multi-sensor  receivers  yields  ambiguity-domain 
direction  finding  and  source  separation  methods  [7]. 

The  purpose  of  this  paper  is  to  consider  stationary,  cyclostationarity  and  nonstationar- 
ity  in  a  unified  framework,  with  applications  to  key  problems  in  array  processing.  In  the 
unified  time-lag  frequency-lag  domain,  cyclostationary  and  nonstationary  signals  demon¬ 
strate  different  properties  which  can  be  used  for  source  selection  and  discrimination,  lead¬ 
ing  to  simplified  approaches  and  improved  performance.  An  example  in  array  processing 
is  provided  where  nonstationary  and  cyclostationary  signals  are  separated  following  simple 
spatial  signature  estimation  exploiting  the  aforementioned  properties. 

II"  Cyclic  Auto-Correlation  Function  and  Ambiguity  Function 

A.  Cyclic  Auto-Correlation  Function 

The  cyclic  auto-correlation  function  pf*  (r)  of  a  scalar  signal  s{t)  is  defined  as 

eiftt5(i  +  r/2)s*(t-r/2)dt  '  (1) 

where  (•)*  denotes  complex  conjugate,  and  r  and  $  denote  the  time-lag  and  frequency-lag 
variables,  respectively.  Signals  with  various  modulations,  such  as  amplitude  modulation 
(AM),  phase  shift  keying  (PSK),  and  frequency  shift  keying  (FSK),  demonstrate  the  cy¬ 
clostationary  property,  that  is,  pfj(r)  /  0  for  non-zero  values  of  $. 
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B.  Ambiguity  Function 

The  ambiguity  function  of  a  signal  s(t)  is  defined  as 

/OO 

ej6ts(t  +  r/2)s*(t  -  r/2)dt  (2) 

-OO 

where  ( 6,t )  are  a  pair  of  frequency-lag  and  time-lag  variables,  respectively.  Comparing  (1) 
and  (2),  it  is  obvious  that  the  cyclic  auto-correlation  function  and  the  ambiguity  function 
are  identical  except  for  a  normalization  factor. 

C.  A  Unified  Framework  for  Discrete-Time  Signals 

In  this  paper,  we  consider  a  unified  framework  for  cyclostationary  and  nonstationary 
signals.  By  considering  N  observations  of  discrete-time  data,  we  define  the  unified  cyclic 
auto-correlation  function  and  (normalized)  ambiguity  function  as 

V.(9,t)  =  fZei«s(t  +  T/2)s-(t-T/ 2).  (3) 

iV  t= 0 

D.  Example  and  Remarks 

In  Fig.  ??,  we  illustrate  the  cyclic  auto-correlation  function  of  a  BPSK  signal,  where 
4  samples  are  generated  from  each  symbol.  128  data  symbols  are  used.  A  rectangular 
waveform  is  used  at  each  symbol  and  the  amplitude  of  the  waveform  is  unity.  For  all 
plots  showing  the  cyclic  auto-correlation  function  (ambiguity  function),  the  time-lag  is 
normalized  by  the  symbol  duration  T,  and  the  frequency  lag  is  normalized  by  the  symbol 
rate  fs  =  1/T.  With  sufficiently  large  number  of  data  samples,  the  cyclic  auto-correlation 
function  has  non-zero  values  only  at  6  =  kfs,  where  k  is  an  integer.  The  cyclic  auto¬ 
correlation  function  depends  on  the  symbol  rate  and  the  pulse  shaping  function. 

In  Fig.  ??,  we  illustrate  the  ambiguity  function  of  a  chirp  signal,  which  has  the  same 
data  length  as  the  BPSK  signal  (note  the  different  r  domain  compared  to  Fig.  1).  The 
start  frequency  of  the  chirp  is  0  and  the  end  frequency  equals  the  symbol  rate  fs  =  1/T 
(the  corresponding  range  of  6  is  from  0  to  n/2).  The  amplitude  of  the  FM  signal  is  1. 
The  ambiguity  function  of  a  chirp  signal  shows  a  linear  signature,  for  which  the  line  of  the 
peak  values  always  passes  through  the  origin  (9  =  r  =  0). 

It  is  clear  that  the  peak  values  of  the  cyclic  auto-correlation  of  a  BPSK  signal  and 
the  ambiguity  function  of  a  chirp  signal  overlap  over  only  a  very  limited  area.  Unlike 
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the  cyclic  auto-correlation  function  of  the  BPSK  signal,  which  has  fixed  and  multiple 
cyclic  frequencies,  the  frequency-lag  for  the  ambiguity  function  of  a  chirp  signal  is  time- 
lag  dependent.  Moreover,  high  values  of  the  ambiguity  function  appear  for  each  time-lag, 
in  contrast  to  the  cyclic  auto-correlation  function  of  a  BPSK  signal  where  the  peaks 
correspond  to  only  specific  time-lag  values.  These  properties  are  very  important  in  source 
discrimination,  as  will  be  discussed  later. 

It  is  noted  that  the  distinctive  signature  contrast  between  the  digital  signal  and  the  FM 
signal  is  not  preserved  in  the  time-frequency  domain,  where  the  autoterm  of  the  digital 
signal  spreads  over  the  entire  time  and  frequency  extent  of  the  signal. 

III.  Array  Processing 

A.  Signal  Model 

When  m  source  signals  impinge  on  an  array  of  n  sensors  (we  assume  that  n>  m),  the 
received  signal  vector  x(t)  =  [xi(t),  •  •  •  ,xn(t)]T  is  expressed  in  the  following  form  ((•)r 
denotes  the  transpose  of  a  matrix  or  a  vector) 

x(i)  =  y  (t)  +  w  (t)  =  As  (t)  +  w  (£)  (4) 

where  y(t)  =  [yi{t),  ■  ■  ■  ,yn{t)]T  is  the  noise-free  signal  vector,  w(f)  =  ■  ■  •  ,wn(t)]T 

is  the  noise  vector,  and  s(t)  =  [si(t),  ■  •  • ,  sm(t)]T  is  the  source  signal  vector.  Moreover, 
A  =  [ai,  •  •  •  ,am]  is  the  n  x  m  mixing  matrix  which  is  unknown.  It  is  assumed  that  the 
noise  vector  w  (t)  is  stationary  and  is  uncorrelated  with  the  source  signal  vector  s  (t)  (the 
noise  vector  may  be  spatially  and  temporally  white  or  colored). 

For  multi-antenna  receivers,  the  ( k ,  l)- th  element  of  the  cyclic  correlation  matrix  (spatial 
ambiguity  matrix)  Ux(0,  t)  can  be  defined  as  the  cyclic  cross-correlation  function  (cross¬ 
ambiguity  function)  between  Xk  and  xi,  expressed  as 

=  '  (5) 

Denote  Us(0,  r)  as  the  cyclic  correlation  matrix  (ambiguity  matrix)  of  the  source  vector  s. 
The  (k,  l)-th  element  of  Us(0,  r)  is  the  cyclic  cross-correlation  function  (cross-ambiguity 
function)  between  Sk(t)  and  si(t)).  We  obtain  the  following  result  from  (5) 

£[Ux(0,r)]  =  AUs(0,r)Aff  (6) 
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for  0/0,  where  (-)H  denote  conjugate  transpose.  Note  that  the  noise  component  is  not 
shown  in  (6)  as  the  expectation  value  of  the  cyclic  auto-  and  cross-correlation  function 
(ambiguity  function)  of  the  stationary  noise  is  zero  for  0  /  0. 

B.  Source  Discrimination  and  Separation 

It  is  noted  that  (6)  is  valid  for  any  point  (0,  r),  0/0.  When  considering  one  BPSK 
and  one  FM  signal,  (6,  r)  regions  can  be  identified  which  include  (0,  r)  points  belonging 
to  only  one  signal.  That  is,  in  the  first  group  Qi  where  the  (0,r)  points  pertain  to  the 
BPSK  signal,  we  have 

E[Ux(0,r)]  =  a1tf,1(0,r)af.  (7) 

Similary,  in  the  second  group  Sl2  where  the  (0,  r)  points  pertain  to  the  FM  signal,  we  have 

JF[Ux(0,r)]  =  a2C/S2(0,r)af.  (8) 

Therefore,  in  either  fix  or  fi2,  £’[UX(0,  r)]  is  rank  one,  and  the  spatial  signature  a,  can 
be  easily  estimated  from  JF[Ux(0,r)]  evaluated  at  any  (0,  r)  €  fl*,,  k  =  1,2.  Given 
finite  period  of  data  observation,  multiple  (0,  r)  points  can  be  incorporated  to  achieve 
an  improved  estimate  of  the  spatial  signature.  A  convenient  way  to  incorporate  multiple 
(0,  r)  points  is  data  augmentation  [3],  which  obtains  the  spatial  signature  as  the  eigenvector 
corresponding  to  the  maximum  eigenvalue  of  VfcVjf,  where  V*  is  an  augmented  spatial 
cyclic  correlation  matrix  (spatial  ambiguity  matrix) 

v,  =  [Ux(0o,  r0),  Ux(0x,  Ti),  •  •  • ,  Ux(0x,_1,  tl—x)],  (9) 

where  Ux(0i,Tj)  G  for  i  =  0, —  1  and  k  =  1,2.  The  high  computations  of 
eigendecomposition  can  be  avoided  by  using  a  simpler,  but  less  noise  robust  approach  in 
which  one,  say,  the  first,  column  of  matrix  VfcVf  is  used  instead. 

Once  the  spatial  signatures  of  both  signals  are  estimated  as  A  =  [ax ,  a2] ,  the  waveform 
of  the  BPSK  signal  can  be  estimated  from 

s  =  A#x  (10) 

where  (■)*  denotes  pseudo-inverse  of  a  matrix. 
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It  is  noted  that,  when  one  of  the  two  signals  is  far  stronger  than  the  other,  the  spatial 
signature  of  the  weak  signal  should  be  estimated  after  the  strong  signal  is  removed  by 
projecting  the  data  to  its  orthogonal  subspace  [8].  An  example  of  recovery  of  a  BPSK 
signal  in  the  presence  of  strong  FM  jammer  is  presented  in  [9]. 


IV.  Simulation  Results 


In  the  simulations,  we  consider  one  BPSK  signal  and  one  chirp  signal  impinging  on 
a  three-sensor  linear  array  with  half-wavelength  spacing.  The  respective  waveforms  are 
described  in  Section  II.  The  input  power  of  both  signals  is  10  dB,  and  the  input  noise 
power  is  0  dB.  The  following  mixing  matrix  is  used: 


A  = 


0.890  +  j0.465  0.962  +  jO.012  ' 

0.065  -  j0.826  0.779  +  jO.481 


-0.761  +  j 0.852  -1.100  -  j0.160 


which  corresponds  to  the  magnitude  value  0.62  of  the  spatial  correlation  coefficient  between 
the  BPSK  and  the  chirp  signal. 

Fig.  ??  shows  the  true  and  estimated  waveforms  of  both  signals.  512  ( 6 ,  r)  ^  0  points 
along  the  chirp  signature  are  used  in  estimating  the  spatial  signature  of  the  chirp  signal, 
whereas  4  (6,  r)  points  at  (±1/T,  ±T/2)  are  used  for  the  BPSK  signal.  The  mean  square 
error  (MSE)  obtained  from  200  independent  trials  is  —15.65  dB  and  —12.51  dB  respectively 
for  the  two  signals,  which  is  very  close  to  the  bound  —15.67  dB  and  —12.66  dB  computed 
using  the  true  mixing  matrix  A  instead  of  A  in  (10).  The  BPSK  signal  has  a  lower  MSE 
of  about  3  dB  because  the  quadrature  noise  component  is  removed. 


V.  Conclusions 

A  generalized  framework  is  presented  for  separation  and  waveform  recovery  of  cyclosta¬ 
tionary  and  nonstationary  signals.  This  approach  is  based  on  the  general  definition  of  the 
cyclic  auto-correlation  function  and  the  ambiguity  function.  Based  on  the  distinct  local¬ 
ization  properties  of  the  nonstationary  signals  and  cyclostationary  signals  in  the  time-lag 
frequency-lag  domain,  source  discrimination  can  be  performed  prior  to  array  processing, 
leading  to  simple  yet  effective  waveform  recovery. 
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Fig.  2.  Time-averaged  ambiguity  function  (magnitude)  of  a  chirp  signal. 
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Fig.  3.  True  ( - )  and  estimated  ( — )  waveforms  of  both  signals. 
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Aperture  Synthesis  for  a  Through-the-wall  Imaging  System 


Fauzia  Ahmad,  Saleem  A.  Kassam,  Gordon  J.  Frazer,  and  Moeness  G.  Amin 


Abstract 

An  aperture  synthesis  scheme  using  subarrays  that  is  based  on  the  coarray  concept  is 
proposed  for  through-the-wall  microwave  imaging  applications.  Simulation  results 
depicting  the  effectiveness  of  the  proposed  synthetic  aperture  technique  for  a  TWI  system 
are  presented.  The  effects  of  incorrect  estimates  of  wall  parameters  and  errors  in  array 
element  placement  on  this  technique  are  also  investigated. 
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1.  Introduction 


The  ability  to  see  through  walls,  doors  and  other  visually  opaque  materials,  using 
microwave  signals  has  become  a  major  area  of  interest  in  a  variety  of  military  and 
commercial  applications.  Through-the-wall  imaging  (TWI)  can  be  used  in  rescue 
missions,  behind-the-wall  target  detection,  surveillance  and  reconnaissance,  and  even 
sensing  through  smoke  and  dust.  Ferris  and  Currie  [1,2]  have  reviewed  existing  and 
under  development  microwave  systems  that  can  detect  the  presence  of  persons  behind 
walls  and  track  their  movement.  Depending  on  the  technology  involved,  these  systems 
can  provide  a  range  resolution  of  a  few  inches  and  are  able  to  sense  moving  targets  and 
measure  their  speed.  However,  most  of  these  systems  have  poor  spatial  resolution. 

Improved  spatial  resolution  can  be  achieved  by  using  an  enlarged  array  aperture. 
However,  with  the  constraints  of  portability  and  low  cost  on  the  system,  an  innovative 
scheme  is  required  for  increasing  the  effective  system  aperture.  In  this  paper,  we  use  an 
aperture  synthesis  scheme  based  on  the  coarray  formalism  for  TWI.  The  concept  of 
coarrays  was  originally  defined  for  narrow-band  far-field  active  imaging  [3,  4],  and  has 
also  been  extended  to  wideband  imaging  [5].  The  coarray  provides  a  convenient  and 
elegant  framework  for  understanding  linear  imaging  techniques.  This  concept  is 
important  for  active  imaging,  as  the  angular  resolution  and  the  point  spread  function, 
which  is  the  response  of  the  system  to  a  point  target,  are  directly  related  to  the  size  of  the 
coarray  of  the  imaging  system. 

The  coarray  completely  characterizes  the  performance  of  an  imaging  system,  and  is 
defined  to  be  the  set  [3] 
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where  St  and  Sr  are  the  sets  containing  the  position  vectors  of  the  elements  in  the 
transmit  and  receive  apertures,  respectively.  It  is  possible  for  two  systems  that  have  the 
same  coarray  to  achieve  the  same  imaging  performance. 

2.  System  Specifications 

In  the  underlying  TWI  system,  we  consider  a  pulsed  radar  system  with  a  bandwidth 
of  666  MHz  centered  at  2  GHz.  The  system  has  a  range  resolution  of  0.225  m. 

2.1  Choice  of  pulse  repetition  frequency  (PRF) 

The  maximum  unambiguous  range  sets  the  upper  limit  on  the  PRF,  whereas  the  lower 
limit  is  set  by  the  desire  to  avoid  Doppler  ambiguities.  For  a  pulsed  radar,  these  limits  on 
PRF /r  are  given  by 

2fd  max  —  fr  —  c/2^«  (2) 

where  Ru  and  /<fmax  are  the  maximum  unambiguous  range  and  maximum  anticipated  target 
Doppler  frequency  respectively.  The  free  space  path  loss  for  the  round  trip  range  of  100m 
at  2  GHz  is  78.4  dB  [6].  In  addition  to  this,  there  will  be  losses  due  to  attenuation  as  a 
result  of  transmission  through  the  wall,  and  further  losses  due  to  the  fact  that  only  part  of 
the  transmitted  signal  gets  reflected  at  each  juncture.  Therefore,  signals  returning  from 
beyond  this  round  trip  range  of  100m  would  be  highly  attenuated.  We  set  the 
unambiguous  range  for  our  TWI  radar  system  to  be  50  m  (one  half  of  the  round  trip 
range),  which  corresponds  to  a  PRF  of  3  MHz.  On  the  other  hand,  considering  the  effect 
of  likely  translation,  rotation,  and  oscillation  of  moving  objects  [7],  we  opt  for  a 


maximum  anticipated  Doppler  frequency  of  250  Hz.  This  sets  the  lower  limit  of  500  Hz 
on  the  PRF.  Accordingly,  500  Hz  <fr  <3  MHz. 


2.2  Effects  of  the  wall  through  which  the  system  is  looking 

The  composition  and  thickness  of  the  wall,  its  dielectric  constant,  and  the  angle  of 
incidence  all  affect  the  strength  and  characteristics  of  the  signal  propagating  through  the 
wall.  The  dielectric  constant  of  a  dry  concrete  wall  lies  between  8  and  12.  The  velocity  of 
a  radio  wave  through  a  medium  of  dielectric  constant  e  is  given  by 


where  c  is  the  speed  of  light.  Therefore,  the  velocity  of  a  signal  propagating  through  a 
concrete  wall  will  decrease  by  a  factor  of  2.8  to  3.4,  inducing  a  bias  in  the  target  range 


measurements  [8], 


3.  Aperture  Synthesis  Using  Subarrays 


In  this  scheme,  we  employ  a  smaller  transmit/receive  array  system  as  a  subarray  to 
synthesize  an  effective  larger  aperture.  The  transmit  and  receive  subarray  pairs  are  used 
to  form  component  complex  images  by  weighted  linear  beamforming,  which  will  then  be 
added  coherently  to  obtain  the  composite  complex  amplitude  image  with  the  desired 
spatial  resolution.  The  composite  image  has  an  effective  coarray  that  is  the  union  of  the 
individual  coarrays  corresponding  to  each  transmit/receive  subarray  pair.  This  and 
several  other  aperture  synthesis  techniques  based  on  the  coarray  have  been  discussed  in 
detail  in  [3]. 
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Since  aperture  synthesis  is  obtained  by  independent  serial  use  of  the  various 
subarrays,  the  overall  system  transmit  and  receive  processing  apparatus  simplifies  to  that 
required  by  one  subarray  combined  with  a  subarray  multiplexer.  The  resolution  of  the 
larger  array  can  thus  be  obtained  by  the  use  of  a  smaller  number  of  processing  channels  at 
the  expense  of  increased  data  acquisition  time. 

3.1  Transmit  and  Receive  Apertures 

Achieving  a  spatial  resolution  of  one-half  a  meter  at  a  range  of  5  meters  along  both 
cross-range  directions  requires  the  length  of  the  coarray  along  both  directions  to  be  at 
least  10X,  where  X=0.15m  is  the  wavelength  corresponding  to  2  GHz.  This  resolution 
could  be  obtained  by  use  of  a  X/2-spaced  11x11  element  square  filled  array  for  both 
transmission  and  reception.  However,  such  an  array  would  require  a  total  of  121 
transmit/receive  elements,  thereby  rendering  it  impractical.  Alternatively,  we  can  design 
transmit  and  receive  array  geometries  of  Fig.  l(a,b).  The  transmit  array  is  a  4-element 
sparse  array  with  an  inter-element  spacing  of  8X.  The  receive  array  is  a  X/2-spaced  8x8 
element  square  array.  The  corresponding  coarray  is  shown  in  Fig.  1(c).  It  has  a  length  of 
11. 5X  along  both  cross-range  directions,  thereby  giving  a  resolution  of  0.43m  at  5m 
range.  In  order  to  apply  the  aperture  synthesis  technique,  we  divide  these  arrays  into 
several  transmit/receive  subarray  pairs.  The  total  number  of  such  pairs  is  determined  by 
the  increased  data  acquisition  time  due  to  serial  operation  of  all  the  subarray  apertures 
and  the  PRF  of  the  system.  In  the  extreme  case  of  a  single  transmitter  and  a  single 
receiver  subarray,  there  will  be  a  total  of  256  subarray  pairs.  This  requires  256  pulses  to 
generate  one  composite  complex  amplitude  image.  If  the  system  is  operating  at  a  PRF  of 
3  MHz,  the  “actual”  PRF  for  the  system  using  subarrays  will  be  3  MHz/256  or  11.72 
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kHz,  which  is  still  23  times  higher  than  the  lower  limit  of  500  Hz  on  the  PRF.  Assuming 
the  scene  is  stationary  over  the  0.80|is  pulse  repetition  interval,  we  opt  for  single 
transmitter,  single  receiver  subarray  apertures  for  the  proposed  system. 

The  subarray  scheme  in  aperture  synthesis  partitions  the  coarray  into  distinct  regions, 
which,  for  our  system,  are  single-element  disjoint  subsets  of  the  whole  coarray.  The 
coarray  subsets  corresponding  to  the  64  subarray  pairs  involving  the  same  transmitter  are 
represented  by  the  same  shade  of  gray  as  that  of  the  transmitter  in  Fig.  1(c). 

3.2  Imaging  System  Model 

We  consider  the  wall  and  the  transmit/receive  subarray  of  omni-directional 
transducers  to  be  located  in  the  x-y  plane,  whereas  the  volume  to  be  imaged  is  located 
along  the  positive  z-axis.  The  volume  of  interest  is  divided  into  a  finite  number  of  pixels 
in  range,  azimuth,  and  elevation.  In  order  to  image  the  set  of  ranges  located  along  a 
particular  elevation  and  azimuth  (0f,cpf),  we  apply  a  time  delay  to  the  transmit  element  to 
focus  it  at  range  Rf  along  (0f,(pf).  Then,  we  transmit  the  pulse,  and  process  the  reflections 
as  follows.  Time  delay  is  applied  to  the  receive  element  so  that  it  is  focused  at  (Rf,0f,(pf). 
The  delayed  output  from  the  receiver  is  matched  filtered  to  detect  I  and  Q  components  of 
the  resulting  signal,  and  a  single  complex  value  I+}Q  is  formed.  This  process  is  repeated 
for  each  range  Rf  in  direction  (0f,<pf)  and  then  repeated  again  for  each  elevation  and 
azimuth  in  the  volume  of  interest.  Component  images  are  formed  by  repeating  this 
process  with  various  subarrays,  and  are  coherently  added  to  form  the  composite  image. 
The  displayed  image  is  the  magnitude  of  the  complex  value  measured  for  each  pixel. 

Let  us  now  consider  in  detail  what  the  value  of  a  single  image  pixel  will  be  with  a 
single  transmitter,  single  receiver  subarray  for  the  case  of  a  single  point  target.  The 
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general  case  of  multiple  targets  can  be  obtained  by  superposition.  Let  the  transmitter  and 
receiver  locations  be  xt  and  xr.  Let  the  target  be  located  in  the  direction 
g=(u,v)=(sin(0)cos((p),  sin(0)sin((p))  at  range  R  as  shown  in  Fig.  2.  Suppose  it  is  desired 
to  image  the  pixel  at  range  Rf  in  the  direction  h=(sin(0f)cos((pf),  sin(0f)sin(cpf)).  As  a  time 
reference,  let  an  element  at  the  origin  be  given  zero  time  delay.  Let 
p(t)=Re { y(t)cxp(jcoct) }  be  the  pulse  transmitted  from  an  element  at  the  origin.  y(t)  is  the 
complex  amplitude  of  the  pulse  at  coc.  Then  the  pulse  received  at  xr  when  the  transmitter 
is  at  xt  is 

q(  t)=Re  { cxp(-a(dt+dr))&(R,g)w,wr  y(t-z)cxp(jcoc(t-x)) }  (4) 

where  a  is  the  attenuation  constant  of  the  wall,  d,  and  dr  are  the  distances  traveled  through 
the  wall  on  transmit  and  receive  respectively,  a(R,g)  is  the  complex  reflectivity  of  the 
target,  w,  and  wr  are  the  weights  applied  to  the  transmit  and  receive  elements  respectively 
and  x  is  the  sum  of  propagation  and  focusing  delays  when  the  array  is  focused  in  direction 
h.  These  delays  incorporate  the  slowing  of  waves  due  to  propagation  through  the  wall. 
The  received  signal  is  passed  through  I  and  Q  matched  filters,  and  their  outputs  are 
sampled  at  time  T  corresponding  to  the  waveform  time  of  flight  for  the  focused  range  Rf. 
Hence, 

I  =  \q{Z)y{S-T)cas{ac{$-T))d$  (5) 

Q  =  -\  q{g)y(g  ~  T)  sin(*>c  (<f  -  T))d%  (6) 

where 

T  =  (2 Rf  ~(d,+dr))/c  +  (dt+dr )/v  (7) 

and  v  is  given  by  (3). 
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3.3  Model  of  the  Walls 


Each  of  the  two  side  walls  and  the  scene  back  wall  is  modeled  as  a  set  of  point 
reflectors  having  the  same  constant  reflectivity  coefficient.  Since  backscattering  from  the 
two  side  walls  is  small  compared  to  the  back  wall,  the  reflectivity  coefficients  as  well  as 
the  density  of  the  point  targets  constituting  the  side  walls  are  given  a  lower  value  relative 
to  the  back  wall.  Figure  3  shows  how  these  walls  would  look  like  in  the  R-u  plane. 

4.  Simulation  Results 

In  this  section,  we  present  simulated  B-Scan  and  C-Scan  images  obtained  using 
subarray  aperture  synthesis  for  the  TWI  system  described  above.  All  of  the  B-scan 
(Range  vs.  u)  images  are  computed  for  v=0.  The  room  being  imaged  is  5m  x  7m.  The 
wall  through  which  the  system  is  looking  is  a  6"  thick  concrete  wall  with  a  dielectric 
constant  s  of  9.  The  room  contains  4  stationary  point  targets,  representing  persons,  whose 
locations  and  reflectivities  are  given  in  Table  1.  Since  the  human  body  is  composed  of 
65%  water,  its  reflectivity  is  quite  high  as  compared  to  that  of  the  walls.  Therefore,  the 
back  wall  is  given  a  reflectivity  of  0.5  whereas  each  side  wall  is  given  a  reflectivity  of 
0.3.  The  reflectivities  are  assumed  to  be  frequency-independent.  Also,  the  one-way 
attenuation  through  the  wall  is  taken  to  be  6  dB  for  6"  thickness  [7,8].  In  all  the  figures, 
we  plot  the  magnitude  of  the  reflectivity  estimate,  with  maximum  value  normalized  to 
unity. 

Figure  4  shows  the  B-scan  image  obtained  under  the  assumption  that  the  wall 
parameters  (thickness  and  dielectric  constant)  are  known  exactly.  We  see  that  our  scheme 
gives  accurate  estimates  of  target  ranges  and  angular  locations  and  is  able  to  resolve  the 
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targets  at  the  same  range.  Figure  5  shows  the  extreme  case  scenario  when  we  assume  to 
have  no  knowledge  of  the  wall  parameters  (estimated  thickness  and  e  are  0  and  1 
respectively).  We  note  that  the  range  estimates  have  large  errors  since  the  slowing  down 
of  the  waves  is  not  accounted  for.  Also  due  to  focusing  errors  in  beamforming,  there  is  an 
error  of  0.1  in  the  u-location  of  the  target  at  3m.  The  u-location  errors  of  the  other  3 
targets  are  not  that  significant.  This  is  because  for  targets  closer  to  the  front  wall  and  at 
angles  way  off  broadside,  the  distance  through  the  wall  is  a  significant  portion  of  the 
target  range  and  hence,  errors  are  more  pronounced.  However,  these  errors  would  occur 
even  if  the  entire  transmit  and  receive  arrays  were  used  for  beamforming  instead  of  the 
subarray  approach.  Due  to  component  image  addition  with  focusing  errors,  the  sidelobes 
are  higher  by  2  dB  approximately  over  those  in  Figure  4. 

Figure  6  shows  the  B-scan  image  obtained  with  exact  knowledge  of  wall  parameters, 
but  with  the  following  errors  in  transmit  array  locations  (2"  in  left,  1 "  in  right,  1 "  in  top, 
0"  in  the  bottom  element).  These  errors  are  possible,  although  large,  in  a  practical  set  up. 
The  beamformer,  however,  assumes  no  such  errors  exist.  Although  the  sidelobes  are 
raised,  no  other  significant  errors  are  observed.  A  C-scan  image  (u  vs.  v)  for  a  fixed  range 
of  6.0375m  under  no  error  conditions  is  shown  in  Figure  7.  Clearly,  the  desired  spatial 
resolution  is  obtained  in  both  cross  range  directions. 

5.  Conclusions 

We  have  examined  the  effectiveness  of  the  aperture  synthesis  technique  using 
subarrays  for  the  through-the-wall  microwave  imaging  problem.  It  is  shown  that  the 
aperture  synthesis  technique  can  be  used  to  achieve  the  desired  cross-range  resolution. 
The  application  of  sparsely  populated  array  geometries  is  desirable  in  TWI  applications 
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to  reduce  system  cost  and  complexity.  We  have  shown  that  the  proposed  scheme  is  robust 
to  a  variety  of  errors  such  as  incorrect  estimates  of  wall  parameters  and  array  element 
location  errors. 
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Figure  1 :  (a)  Transmit  array;  (b)  Receive  array;  (c)  Coarray 
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Figure  2:  Path  of  a  signal  transmitted  from  xt,  reflected  from  xs,  and  received  at  xr.  The 
wall  width  is  along  x-axis  and  y-axis  represents  its  height,  with  the  origin  at  the 
center  of  the  wall.  The  length  of  the  room  lies  along  the  positive  z-axis. 


Figure  3:  (a)  Geometry  depicting  the  ranges  and  angles  corresponding  to  points  on  the 
side  and  the  opposing  walls  in  the  x-z  plane;  (b)  The  side  and  the  opposing  wall 
in  the  R-u  plane. 
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Figure  7:  C-Scan  image  at  range  6.0375  m. 


68 


IEEE  Sensor  Array  and  Multichannel  Signal  Processing  Workshop, 

Rosslyn,  VA,  August  2002. 


Time-frequency  ESPRIT  for 
direction-of-arrival 
estimation  of  chirp  signals 


Aboulnasr  Hassanien,  Alex  B.  Gershman,  and  Moeness  G.  Amin 


Abstract 

In  this  paper,  we  develop  an  ESPRIT-type  algorithm  for  estimating  the  Directions-Of- Arrival  (DOA’s) 
of  multiple  chirp  signals  using  Spatial  Time-Frequency  Distributions  (STFD’s).  An  averaged  STFD  matrix 
(or  multiple  averaged  STFD  matrices)  are  used  instead  of  the  covariance  matrix  to  estimate  the  signal  and 
noise  subspaces.  The  proposed  algorithm  is  shown  to  provide  a  significant  performance  improvement  over 
the  traditional  ESPRIT  algorithm  for  FM  sources,  specifically  in  situations  with  closely  spaced  sources 
and  low  Signal-to-Noise  Ratios  (SNR’s).  Simulation  results  are  provided  to  illustrate  the  performance  of 
the  proposed  approach  in  scenarios  with  multiple  narrowband  chirp  sources. 
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I.  Introduction 


The  problem  of  DOA  estimation  in  the  presence  of  chirp  and  FM  signal  sources  in  sensor 
arrays  arises  in  Synthetic  Aperture  Radar  (SAR),  Synthetic  Aperture  Sonar  (SAS),  inverse 
SAR  and  SAS  (ISAR  and  ISAS),  Doppler  radar  and  sonar  imaging,  as  well  as  in  mobile 
communications,  where  FM  signal  waveforms  can  be  intentionally  transmitted  [l]-[9].  In 
conventional  array  processing,  subspace  methods  (for  example,  MUSIC  and  ESPRIT)  are 
commonly  applied  and  achieve  excellent  performance  at  a  moderate  computational  cost. 
However,  conventional  subspace  methods  are  based  on  the  assumption  that  the  received 
signals  are  stationary.  As  a  result,  the  performance  of  conventional  subspace-based  DOA 
estimation  techniques  can  degrade  when  dealing  with  nonstationary  chirp  signals. 

Recently,  STFD’s  have  emerged  as  an  efficient  means  of  array  processing  in  the  case  of 
multiple  FM  signals  [10]-[14].  Spreading  the  noise  power  while  localizing  the  sources  in  the 
time-frequency  domain  helps  to  improve  the  DOA  estimation  performance  and  enhance 
robustness  against  sensor  noise.  However,  the  existing  narrowband  STFD-based  DOA 
estimation  techniques  [10]-[13]  are  based  on  a  spectral  search  and,  as  a  result,  have  high 
computational  costs. 

In  this  paper,  a  search-free  STFD-based  (time-frequency)  ESPRIT  algorithm  is  devel¬ 
oped.  In  order  to  obtain  improved  estimates  of  the  signal  and  noise  subspaces,  we  use 
an  averaged  STFD  matrix  (or  multiple  averaged  STFD  matrices)  in  place  of  the  covari¬ 
ance  matrix  which  is  used  in  the  conventional  ESPRIT  algorithm.  This  averaging  involves 
time-frequency  points  that  correspond  to  source  signatures  with  a  maximal  energy.  The 
source  DOA’s  are  then  estimated  using  either  the  least  squares  (LS)  or  the  total  least 
squares  (TLS)  ESPRIT  approach  [16]- [18]. 

The  proposed  technique  enables  to  separate  the  signals  in  different  averaged  STFD 
matrices  prior  to  DOA  estimation  and,  therefore,  makes  it  possible  to  estimate  the  source 
DOA’s  in  the  case  when  the  number  of  array  sensors  is  less  than  the  number  of  sources. 
Moreover,  closely  spaced  sources  with  well  separated  time-frequency  signatures  can  be 
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efficiently  resolved  by  separating  them  in  different  averaged  STFD  matrices  and  applying 
the  ESPRIT  algorithm  to  each  of  them  independently. 

To  validate  the  effectiveness  of  the  proposed  technique  and  compare  its  performance 
with  the  conventional  array  processing  methods,  simulation  results  for  chirp  sources  are 
presented.  It  is  shown  that  significant  performance  gains  can  be  achieved  by  the  proposed 
algorithm  compared  to  the  conventional  LS-  and  TLS-based  ESPRIT  techniques. 

II.  Signal  Model 

Consider  L  chirp  signals  impinging  on  a  Uniform  Linear  Array  (ULA)  of  M  sensors. 
The  array  output  (data  snapshot)  vector  x(t)  G  (CMxl  is  modeled  as 

x(t)  =  Ad(t)  +  n(t)  ,  t  =  l,...,N  (1) 

where  d (t)  €  CLxl  and  n(t)  G  <CMxl  are  the  vectors  of  signal  waveforms  and  sensor  noise, 
respectively,  N  is  the  number  of  data  snapshots  available,  and  A  G  CMxi  is  the  array 
direction  matrix.  It  is  assumed  that  A  is  full  rank  (which  means  that  the  steering  vectors 
corresponding  to  L  different  angles  of  arrival  are  linearly  independent).  We  assume  that 
the  noise  is  a  spatially  and  temporally  white  zero-mean  process,  i.e. 

E{n(*)n«(s)}=<72I,5M  (2) 

where  I  is  the  identity  matrix,  a2  is  the  noise  variance,  St<s  is  the  Kronecker  delta,  and 
(■)H  stands  for  the  Hermitian  transpose.  The  direction  matrix  can  be  expressed  as 

A=[a(01),...,a(0i)]  (3) 

where 

a(0)  =  [1,  ej^Asin9, . . . ,  ^^M~l'>sind]T  (4) 

is  the  steering  vector,  {Ok}k= i  are  the  signal  DOA’s,  uj  is  the  central  frequency,  c  is  the 
propagation  speed,  A  is  the  interelement  spacing,  and  (-)T  denotes  the  transpose.  Note 
that  this  model  corresponds  to  the  assumption  of  narrowband  chirp  signals  where  changes 


of  the  central  frequency  within  the  observation  interval  are  negligibly  small  (i.e.,  the  matrix 
A  is  time  independent)  [10]-[13].  The  case  of  wideband  chirp  signals  is  addressed  in  [14] 
and  [15]. 

The  sample  Spatial  Pseudo-Wigner-Ville  Distribution  (SPWVD)  matrix  is  given  by 
[10]- [13] 

(*--l)/2 

Dxx(i,  /)  =  E  x(*  +  r)xH(*  “  r)e~3i7rfT  (5) 

t=-(K- l)/2 

where  K  is  the  odd  window  length.  Taking  the  expectation,  and  assuming  that  the  source 
waveforms  and  sensor  noise  are  statistically  independent,  we  have  that  [10],  [14] 

Dxx(t,/)  =  E{Dxx(t,/)} 

=  ADdd(t,/)AH  +  (72I  (6) 

where 

Ddd(f, /)  =  E{Ddd  (£,/)}  (7) 

(K- 1)/2 

Ddd(t,/)=  E  d(t  +  r)dH(t-r)e-^  (8) 

r=-(K- 1)/2 

The  relationship  (6)  holds  true  for  each  ( t ,  /)  point.  In  order  to  reduce  the  effect  of  the 
sensor  noise,  an  averaging  over  multiple  ( t ,  /)  points  (corresponding  to  the  autoterms  of 
the  STFD)  can  be  used.  The  eigenstructure  properties  of  the  averaged  STFD  matrix 
can  be  exploited  to  estimate  the  signal  DOA’s  in  a  similar  way  as  in  the  conventional 
subspace- based  array  processing  techniques  [10]- [14]. 

In  practical  situations,  the  sample  STFD  matrices  (5)  are  used  instead  of  the  exact 
(statistically  expected)  matrices  (6).  In  the  case  of  sources  with  distinct  time-frequency 
signatures,  it  is  possible  to  construct  the  averaged  STFD  matrices  over  time-frequency 
points  belonging  to  a  subset  of  the  sources  in  the  field-of-view.  Using  this  approach, 
sources  with  close  angular  spacing  can  be  efficiently  resolved  by  separating  them  in  dif¬ 
ferent  averaged  STFD  matrices  and  applying  a  DOA  estimation  algorithm  to  each  matrix 
independently. 
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III.  Time-Frequency  ESPRIT 


Let  the  averaged  STFD  matrix  Dxx  be  formed  by  averaging  of  multiple  SPWVD  matri¬ 
ces  computed  at  PJ  different  time-frequency  points  { tp ,  /*}  (p  =  1, . . . ,  P;  i  =  1, . . . ,  J) 
that  belong  to  the  time-frequency  signatures  of  L0  source  signals  (L0  <  L  <  M): 

£>x*  =  ££d*x(W<)  (9) 

P~1  2=1 


Note  that  the  values  of  P  and  J  may  vary  with  time  and  can  be  determined  by  means  of 
detection  of  the  source  time-frequency  signatures,  see  [14],  [19],  and  [20].  Let  us  divide  a 
ULA  of  M  sensors  into  two  identical  subarrays  of  M  —  1  sensors  shifted  by  the  interelement 
spacing  A,  as  shown  in  Fig.  1.  As  in  the  conventional  ESPRIT  case,  define  the  sub¬ 
matrices  Ai  and  A2  by  deleting  the  last  and  first  rows  from  A  respectively,  i.e.  let 


A  = 


Then,  Ai  and  A2  are  related  as 


Ai 

first  row 

last  row 

A2 

A  2  =  Ax$ 


(10) 


(11) 


where 

$  =  diagle7*'1, . . . ,  e^o} 


(12) 


and 


Hi  =  (a>/c)Asin0, 


(13) 


are  the  source  spatial  frequencies. 

Let  Us  be  the  matrix  formed  from  the  eigenvectors  of  Dxx  that  correspond  to  the  L0 
eigenvalues  with  the  largest  magnitude.  As  the  columns  of  the  steering  matrix  A  and  the 
matrix  Us  span  approximately1  the  same  (signal)  subspace  of  the  dimension  L0,  there 
exists  a  nonsingular  L0  x  L0  matrix  T  such  that 


Us  AT 


'Note  that  Us  is  obtained  from  the  eigendecomposition  of  the  sample  STFD  matrix  Dxx. 


(14) 
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Applying  this  transformation  to  the  sub-matrices  Ai  and  A2,  we  obtain 

U5,i  *  AiT 
Uc2  -  A2T 


(15) 

(16) 


where  U54  and  U5)2  are  formed  by  deleting  the  last  and  first  rows  of  Us,  respectively, 
that  is 

Us  =  |  |  =  |  |  (17) 


Us.i 

first  row 

last  row 

Us, 2 

Using  (11)  and  (15)-(16)  yields 

Us, 2  —  Us.i^ 

where  the  matrices  9/  and  are  related  as 

This  means  that  are  the  eigenvalues  of  the  matrix  \I>  [16]-[18]. 

Now,  we  can  formulate  our  time-frequency  ESPRIT  algorithm  as  follows: 


(18) 


(19) 


•  Step  1:  Compute  the  sample  SPWVD  matrices  Dxx(t,  /)  for  all 
time-frequency  points  of  interest  and  select  the  maximal  energy 
points  that  belong  to  the  source  signatures. 

•  Step  2:  Compute  the  averaged  STFD  matrix  Dxx  for  a  previously 
selected  part  of  sources  (for  Lq  <L  sources)  by  means  of  involving 
in  the  averaging  process  only  the  time-frequency  points  that  belong 
to  the  spatial  signatures  of  these  L0  sources. 

•  Step  3:  Compute  the  eigendecomposition  of  Dxx  and  obtain  Us- 

•  Step  4:  Obtain  an  estimate  of  the  matrix  \E'  by  solving  (18) 
using  either  LS  or  TLS  approaches. 

•  Step  5:  Obtain  estimates  of  the  spatial  frequencies  from  the 
eigenvalues  of  4r  and  use  them  to  find  the  estimates  of  the  source 
DOA’s  9i. 

•  Step  6:  Repeat  Steps  2-5  for  other  (selected)  parts  of  sources. 
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Note  that  algorithms  are  available  to  classify  auto-  and  cross-terms  in  STFD’s  [19],  [20]. 
These  techniques  can  be  used  in  Step  2. 

The  possibility  of  separating  sources  in  different  averaged  STFD  matrices  prior  to  DOA 
estimation  can  essentially  improve  the  performance  of  the  ESPRIT  algorithm,  especially 
in  the  low  SNR  case  and  in  the  presence  of  closely  spaced  sources  which  are  well  separated 
in  the  time-frequency  domain.  However,  this  will  increase  the  computational  costs  relative 
to  the  conventional  ESPRIT  algorithm  because  in  this  case  ESPRIT  should  be  applied 
simultaneously  to  several  averaged  STFD  matrices. 

IV.  Simulation  Results 

We  assume  a  ULA  with  A  =  A/2.  This  array  is  divided  into  2  subarrays  with  the  inter¬ 
subarray  displacement  A/2  as  shown  in  Fig.  1.  Two  narrowband  chirp  signals  impinge 
on  the  array  from  the  sources  located  at  6\  =  3°  and  02  =  6°.  After  downconversion,  the 
source  waveforms  are  modeled  as 

di(t)  =  eJ(«>it+0it2/2) 

d2(t)  =  e^+/32t2/2) 

The  initial  discrete-time  frequencies  of  the  source  signals  are  chosen  to  be  u>\  =  0  and 
u>2  =  7T  while  their  chirp  rates  are  assumed  to  be  =  0.002  and  /?2  =  —0.002,  respectively. 
An  observation  interval  of  N  =  255  snapshots  is  considered.  Figure  2  shows  the  pseudo- 
Wigner-Ville  distribution  of  the  signals  in  the  first  array  sensor.  The  noise  is  modeled 
as  a  complex  Gaussian  zero-mean  spatially  and  temporally  white  process.  A  total  of  300 
independent  Monte-Carlo  simulation  runs  have  been  used  to  obtain  each  simulated  point. 
The  averaged  STFD  matrix  Dxx  is  computed  for  each  source  signal  separately  by  averaging 
the  sample  STFD  matrices  computed  at  150  different  time-frequency  points  that  belong 
to  the  source  signatures. 

In  the  first  example,  we  use  an  array  of  M  =  10  sensors.  Figure  3  displays  the  DOA 
estimation  Root-Mean-Square  Errors  (RMSE’s)  versus  the  SNR  for  the  conventional  LS- 
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and  TLS-ESPRIT  algorithms,  as  well  as  for  the  proposed  time-frequency  modification 
of  the  LS-  and  TLS-ESPRIT  techniques.  Additionally,  the  so-called  deterministic  CRB 
[21]  is  shown  in  this  figure  (the  latter  bound  is  computed  under  the  assumption  that  the 
source  waveforms  are  unknown  deterministic  sequences).  ^From  Fig.  3,  it  is  clear  that 
the  time-frequency  ESPRIT  algorithm  has  a  substantial  improvement  over  conventional 
ESPRIT.  This  is  especially  true  in  the  low  SNR  case.  Note  that  neither  time-frequency 
ESPRIT  nor  conventional  ESPRIT  approaches  the  deterministic  CRB. 

In  the  second  example,  the  source  SNR  is  fixed  and  equal  to  4  dB,  while  the  number 
of  array  sensors  is  varied.  Figure  4  shows  the  RMSE’s  of  the  same  four  methods  as  in 
the  previous  example  versus  the  number  of  sensors  M.  It  is  evident  from  this  figure 
that  the  time-frequency  ESPRIT  algorithm  has  substantially  better  performance  than  the 
convectional  ESPRIT  technique.  Furthermore,  at  small  values  of  M,  the  RMSE’s  of  time- 
frequency  ESPRIT  can  be  lower  than  the  deterministic  CRB.  To  explain  this  behavior,  we 
note  that  the  displayed  bound  assumes  unknown  deterministic  signal  waveforms  [21]  and, 
therefore,  does  not  take  into  account  the  chirp  signal  structure  [8].  However,  this  structure 
is  exploited  in  the  proposed  algorithm.  Obviously,  this  can  cause  situations  when  the  DOA 
estimation  accuracy  of  the  proposed  technique  is  better  than  the  deterministic  CRB. 

V.  Conclusions 

A  time-frequency  ESPRIT  algorithm  is  introduced  for  DOA  estimation  of  chirp  signals 
in  sensor  arrays.  The  proposed  algorithm  is  based  on  the  concept  of  Spatial  Time-Frequen¬ 
cy  Distributions  (STFD’s)  and  employs  multiple  averaged  STFD  matrices,  instead  of  the 
covariance  matrix  (used  in  conventional  array  processing  methods),  to  obtain  the  estimates 
of  the  signal  DOA’s.  Computer  simulations  show  that  in  scenarios  with  chirp  signals,  the 
proposed  technique  outperforms  the  conventional  ESPRIT  algorithm.  The  performance 
improvement  is  especially  high  in  the  case  when  the  SNR  is  low  or  the  sources  are  closely 
spaced. 
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Fig.  1.  The  array  structure. 


Fig.  2.  Pseudo- Wigner-Ville  distribution  of  the  source  waveforms. 
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Fig.  3.  The  DOA  estimation  RMSE’s  versus  the  SNR.  First  example:  dx  -  3°,  02  =  6°,  and  M  =  10. 
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Fig.  4.  The  DOA  estimation  RMSE’s  versus  the  number  of  sensors  M.  Second  example:  6\  =  3°,  62  =  & 
and  SNR  =  4  dB. 
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Abstract 

The  problem  of  selecting  auto-  and  cross-terms  of  time-frequency  distributions  (TFDs)  of  nonstation¬ 
ary  signals  impinging  on  a  multi- antenna  receiver  is  considered.  A  detection  approach  is  introduced  which 
allows  performance  measurement  and  comparison  of  various  schemes  via  receiver  operating  characteristics. 
Array  averaging  and  array  differencing  techniques  are  both  employed  to  form  a  basis  for  time-frequency 
(t-f)  point  selection.  The  proposed  classification  method  is  evaluated  against  the  bootsrap-based  method. 
It  is  shown  that  the  former  offers  improved  performance  and  simplified  implementations. 
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I.  Introduction 


Time-frequency  distributions  have  recently  been  proposed  for  applications  to  array  sig¬ 
nal  processing  problems  [1],  [2].  For  this  purpose,  spatial  time-frequency  distributions 
(STFDs)  have  been  introduced  and  represented  in  a  matrix  form.  The  elements  of  a  STFD 
matrix  are  the  time-frequency  distributions  and  the  cross-time  frequency  distributions  of 
the  data  received  at  the  multi-sensor  array.  It  has  been  shown  that  the  relationship  be¬ 
tween  the  spatial  time-frequency  distributions  of  the  sensor  data  and  the  time-frequency 
distributions  of  the  sources  is  identical  to  that  of  the  sensor  data  covariance  matrix  and  the 
sources’  correlation  matrix.  This  key  property  permits  direction  finding  and  blind  source 
separation  to  be  performed  using  the  sources’  time-frequency  localization  properties. 

Blind  source  separation  has  been  typically  solved  using  statistical  information  available 
on  source  signals,  including  second  or  higher  order  statistics.  A  primary  contribution  in 
this  area  was  to  show  that  the  STFDs  are  an  effective  alternative  to  separating  sources 
whose  signatures  are  different  in  the  t-f  domain  [1],  [2].  Successful  applications  of  STFDs 
to  source  separation  require  computing  STFDs  at  different  t-f  points.  The  auto-terms, 
which  enforce  the  diagonal  structure  of  the  source  TFD,  are  then  incorporated  into  a  joint- 
diagonalization  (JD)  technique  to  estimate  the  mixing,  or  the  array  manifold  matrix.  This 
technique  was  further  generalized  to  incorporate  cross-term  STFD  matrices  by  performing 
combined  JD  and  joint  off-diagonalization  (JOD)  [3].  In  both  of  these  methods,  proper 
selection  of  either  auto-  or  cross-term  locations  in  the  time-frequency  plane  has  remained 
an  important  and  necessary  task. 

A  statistical  test  to  decide  whether  a  t-f  location  is  dominated  by  auto-  or  cross-terms 
has  been  proposed  in  [4].  Array  averaging  of  the  TFDs  is  used  to  reduce  the  cross-terms 
without  smearing  the  auto- terms  [5],  allowing  the  auto-terms  to  be  more  pronounced  and 
easier  to  detect  in  the  t-f  domain.  However,  as  shown  in  [6],  there  are  advantages  to  using 
both  auto-  and  cross-terms  in  combined  JD/JOD.  This  requires  devising  an  equivalent 
technique  for  cross-term  enhancement  and  selections. 
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In  this  paper  we  propose  a  simple  method  for  the  suppression  of  auto-terms  and  en¬ 
hancement  of  cross-terms.  A  procedure  for  automatic  selection  of  auto-  and  cross-  terms 
is  proposed  and  a  performance  comparison  to  the  bootstrap  based  scheme  is  given. 

II.  Signal  Model 

We  consider  the  general  problem  of  a  number  of  independent  non-stationary  sources 
impinging  on  an  array  of  sensors  which  has  at  least  as  many  elements  as  the  number  of 
sources.  The  model  for  the  array  output  is  given  as  follows, 

x(t)  =  As(t)  +  n(t),  (1) 

where  x(t)  =  [x^t), . . .  ,  xm(t)]'  is  the  array  output  vector,  s(t)  =  [sx (t), . . .  ,  sd(t)]'  is  the 
dxl  source  signal  vector  and  n(t)  is  the  mx  1  independent  and  identically  distributed  noise 
vector  of  variance  a2,  at  time  t.  The  mxd  mixing  matrix  A  is  of  the  form  A  =  [ax, . . .  ,  a<i], 
where  a;  =  [an, . . .  ,  ami]'  denotes  the  steering  vector  for  source  i.  Here  m  represents  the 
number  of  array  elements  and  d  denotes  the  number  of  source  signals. 

A.  Spatial  TFD 

The  discrete  form  of  Cohen’s  class  of  TFD’s,  for  a  signal  x(t)  is  given  by  [7] 

OO 

Dxx(t ,  /)  =  l)x(t  +  m  +  l)x*(t  +  m- 

Z,m=— oo 

where  0(m,  /)  is  the  kernel  defining  the  distribution.  The  STFD  matrix,  based  on  Dxx(t,  f ) 
has  been  introduced  in  [1]  and  is  given  by 

OO 

Dxx(t,  f)  =  l)x(t  +  m  +  l)xH(t  +  m  -  l)e~j4,rfl. 

l,m=— oo 

In  the  (noise-free)  case  of  auto-sensor  TFD’s, 

DxkXk  —  X^u=l  Sv=l  akuakvDsusv(t)  f)  (2) 

where  H  is  the  Hermition  operator  and  DSiSj (t,  /)  is  the  ij  element  of  Dss(t,  f);  the  source 
TFD  matrix.  The  ij  element  of  matrix  Dxx(t,  f)  is  equal  to  DXiXj(t ,  /),  which  is  obtained 
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from  the  expression  for  Dxx(t,f )  by  replacing  the  bilinear  product  of  x(t)  with  the  data 
received  at  sensors  i  and  j. 

Below  we  consider  only  the  pseudo  Wigner-Ville  distribution  [7]  with  odd  window  length 
L  and  henceforth  refer  to  the  STFD  matrix  as  Dxx  where  the  (t,  f )  parameters  are  dropped 
for  convenience.  The  same  applies  to  Dss. 

B.  Test  Statistic 

In  [4],  identification  of  auto-  and  cross-terms  is  achieved  using  the  test  statistic 

trace  {D^}  j  norm  ,  (3) 

where  trace  {•}  is  the  matrix  trace  and  norm  {•}  is  the  Frobenius  norm  of  a  matrix.  With 
a 2  being  an  estimate  of  a2, 

Dzz  =  W(DXX  -  d2I)WH,  (4) 

where  W  is  an  estimate  of  the  whitening  transform  W  such  that  WADSSAHWH  is  a 
unitary  transformation  of  Dss. 

It  is  proposed  to  discriminate  between  noise  and  either  auto-  or  cross-terms  by  using  the 
variance  of  the  test  statistic.  As  only  a  single  value  of  the  test  statistic  is  known  at  any 
t-f  point,  the  variance  is  estimated  using  a  bootstrap  resampling  scheme  by  [4].  Once  the 
noise  regions  in  the  t-f  domain  are  identified,  a  threshold  is  set  to  discriminate  between 
the  auto-  and  cross-term  locations. 

III.  Selection  of  Auto-terms 

We  denote  the  set  of  all  t-f  locations  which  contain  the  auto-terms,  by  Sa  ■  The  test  is 
formulated  as  follows: 

Ho : 

Ha : 


{to,fo)  £  Sa  , 


(5) 
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where  the  null  hypothesis  Ho  is  that  the  point  (to,  /o)  is  not  an  auto-  term.  To  apply  this 
test,  we  define 


j(t,/)  G  ST  :  ^  trace (Dss(t,f)}  >  <yA J 


(6) 


where  St  is  the  set  of  all  locations  in  the  t-f  plane  and  7^  is  an  arbitrarily  selected  bound. 


A.  Array  Averaged  TFD 


The  array  averaged  TFD  improves  signal  synthesis  performance  by  suppressing  the 
cross-terms  [5].  It  is  given  by 


D(t,f) 


m  trace  {Dxx} 
Eli  Eli  PuvDs. 


(<>/)> 


(7) 


where  fduv  =  ^au  is  the  spatial  correlation  coefficient  of  sources  u  and  v.  The  value 
of  /3UV  is  one  for  auto-source  terms  and  less  than  or  equal  to  one  for  cross-source  terms. 
This  coefficient  takes  a  zero  value  if  the  respective  sources  u  and  v  have  orthogonal  spatial 
signatures. 

The  whitened  array  average  is  defined  as 


Dz  =  trace  {WDXXWH}  /d  (8) 

with  expected  value  E[  Dz]  =  ( trace  {Dss}  +  a 2  7w  )  /d,  where  7w  =  trace  {WWH}. 
This  indicates  that  a  test  statistic  useful  for  detection  of  auto-terms  can  be  based  on  the 
whitened  array  average  since  total  cancellation  of  cross-terms  is  achieved. 


B.  Setting  the  Threshold 

In  the  case  of  FM  signals  which  have  constant  power  over  the  observation  interval,  the 
variance  of  the  whitened  array  average  can  be  approximated  by 

Var[  Dz]  =  »  (E)'  L  (2  Tw  +  Tw*) , 
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where  7w2  =  trace  {(WWH)2}.  As  Dz  is  approximately  normally  distributed,  a 
threshold  7  for  classification  of  auto-terms  is  set  according  to 

7  =  7 A  +  v2  Tw/d  +  Nx-aO  oz, 


where  Na  is  the  inverse  of  the  standard  normal  distribution  at  probability  a. 


IV.  Selection  of  cross-terms 

We  denote  the  set  of  all  t-f  locations  that  contain  cross-terms  of  interest,  by  Sc  •  The 
test  is  formulated  as  follows: 


Ho '  (to,  /o)  ^  Sc 

Ha  ■  (to,  /o)  €  Sc 
where  the  set  Sc  is  defined  in  a  similar  manner  to  (6). 


(9) 


A.  Array  Differenced  TFD 

In  order  to  suppress  the  auto-terms  while  supporting  the  cross-terms,  we  define  the 
array  differenced  TFD 

-  m 

b(t,  f )  =  2Z  (0«.  («,/)-  D(t,  f)f  (10) 

u=l 

which  corresponds  to  the  sample  variance  of  the  auto-sensor  TFD’s.  Expanding  the  sum¬ 
mand  with  (2)  and  (7)  we  obtain 


f)  —  DXkXk(t,  f)  D(t,  f) 


d  d  /  ^  m  \ 

= y ( akuakv  ~  ~  y  ^  ^lu^iv }  dSuSv  ( t ,  /) 

u=l  V  m  1=1  / 


and  for  u  =  v  (auto-terms) 


Here  auv  corresponds  to  the  complex  weighting  of  source  v  at  element  u  of  the  array.  The 
array  difference  results  in  complete  cancellation  of  auto-terms,  provided  each  sensor  of  the 
array  has  the  same  gain,  even  if  the  sources  have  different  power.  The  array  difference 
method  is  illustrated  in  Figure  VII,  where  there  are  two  chirp  signals  present  having  DOA’s 
of  10  and  —10  degrees  at  a  2  sensor  array. 

B.  Setting  the  Threshold 

The  variance  of  the  array  difference  is  signal  dependent  and  not  easily  approximated 
unless  the  SNR  is  very  low.  A  threshold  to  automatically  detect  the  cross-terms  may  be 
set  by  the  following  function 

7  =  /i  +  /?<7  (11) 

where  fi  and  a  are  estimates  of  the  mean  and  standard  deviation  of  the  array  difference 
respectively.  Due  to  the  peaks  of  auto-  or  cross-terms  in  the  distribution,  the  mean  may  be 
estimated  using  the  median  of  each  time-slice,  and  the  standard  deviation  as  the  average 
absolute  deviation  from  this  value.  A  suitable  value  for  the  parameter  /3  may  be  in  the 
range  3  to  5  depending  on  the  rate  of  false  classification  that  is  considered  acceptable. 

V.  Point  Classification 

In  blind  source  separation  (BSS)  and  direction  of  arrival  (DOA)  estimation,  the  ultimate 
goal  of  point  selection  is  to  allow  construction  of  STFD  matrices  with  either  a  strong 
diagonal  or  off-diagonal  structure.  This  implies  that  only  peaks  of  the  auto-  and  cross- 
TFD’s  should  be  in  Sa  or  Sc  ■  In  the  case  of  a  signal  synthesis  problem,  we  wish  to 
select  as  many  of  the  desired  signal’s  auto-TFD  terms  as  possible  without  including  cross- 
or  noise-terms.  This  will  determine  the  appropriate  value  for  the  bounds  7 a  or  7 c- 

A.  Point  Classification  Scheme 

There  are  t-f  locations  dominated  by  either  auto-,  cross-  or  noise-terms.  Using  the 
array  averaged  (AA)  TFD  cross-terms  are  suppressed  while  the  array  differenced  (AD) 
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TFD  suppresses  auto-terms.  This  translates  the  classification  problem  with  three  classes, 
into  two  classification  problems  each  containing  two  classes.  A  proposed  point  classification 
scheme  is  illustrated  in  Figure  2. 

The  method  of  point  classification  involves  thresholding  of  the  appropriate  TFD  to  ob¬ 
tain  a  set  of  t-f  points.  The  TFD  is  also  smoothed  using  a  two-dimensional  filter  to  reduce 
the  effect  of  noise,  and  a  threshold  is  applied  to  obtain  another  set  of  t-f  points.  This  set  of 
points  is  processed  using  binary  morphological  operations  under  the  assumption  that  auto- 
and  cross-terms  will  have  some  structure  or  connected  shape  in  the  TF  plane.  The  block 
labelled  O-C  in  Figure  2  signifies  opening  and  closing  operations  [8]  are  used.  The  final 
decision  is  obtained  as  the  intersection  of  these  sets  of  t-f  points.  The  result  is  a  reduction 
in  the  false  classification  rate  compared  to  that  observed  with  simple  thresholding. 

B.  Performance  Evaluation 

Performance  of  the  classification  scheme  is  evaluated  by  means  of  an  operating  charac¬ 
teristic  (OC)  curve.  This  gives  the  probability  of  (correct)  classification  for  a  particular 
probability  of  false  classification.  The  set  of  t-f  locations  obtained  by  point  classification 
of  auto-  and  cross-terms  are  respectively  denoted  Sa  and  <SC  .  Accordingly,  empirical 
probabilities  of  detection  and  false  alarm: 

Pd  =  #{  $4  p|  Sa  }  /  #{  SA  }  , 

Pfa=  #{  SA  -  SA  }/#{&•-  SA  }  , 

where  #{•}  denotes  the  number  of  elements  in  the  set. 

VI.  Simulation 

In  the  following  simulations  we  choose  7 4  and  7 c  to  be  3dB  below  the  peak  of  the  source 
signals  auto-  and  cross-terms  respectively.  The  threshold  of  each  point  detection  scheme 
is  varied  in  order  to  generate  the  operating  characteristic.  Herein  point  detection  using 
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the  array  average  and  array  difference  will  be  referred  to  as  method  one  and  the  bootstrap 
resampling  scheme  will  be  called  method  two. 

In  the  first  case,  we  consider  a  scenario  with  two  source  signals  impinging  on  a  three 
sensor  array  at  0  dB  SNR.  The  source  signals  Si(t)  and  s2(t)  are  chirps  with  start  and  end 
frequencies  (0.37r,0.l7r)  and  (0.257T,  O.OItt)  respectively.  The  OC’s  of  method  one  and  two 
are  shown  in  Figure  3. 

In  the  second  case  we  consider  a  scenario  with  four  source  signals  impinging  on  a  five 
sensor  array  at  5  dB  SNR.  The  source  signals  Si(t)  and  s2(t)  are  as  defined  previously. 
Signal  s3(t)  is  a  chirp  with  start  and  end  frequency  of  (0.27T,  0.3257t).  The  signal  s4(t) 
is  of  the  form 


s4(t)  =  exp[j27r(Kit1,5  +  K2t)], 

where  K\  and  K2  are  chosen  such  that  the  start  and  end  frequencies  are  (0.257T,  0.57t). 
The  signal  s4(t)  has  a  non-linear  instantaneous  frequency  and  therefore  the  auto-terms 
are  not  properly  localized,  when  using  WVD.  The  OC’s  of  method  one  and  two  are  shown 
in  Figure  4. 

A.  Discussion 

Simulations  show  that  in  low  SNR,  correct  identification  of  both  auto-  and  cross-terms 
is  achievable  with  probability  of  false  classification  less  than  2%.  In  the  first  scenario,  only 
30%  of  the  total  cross-terms  were  identified  for  2%  miss-classification.  However,  for  the 
application  of  BSS,  only  a  few  strong  cross-terms  locations  are  required. 

Similar  performance  in  terms  of  operating  characteristic  has  been  achieved  using  both 
methods  for  point  classification.  In  terms  of  computational  complexity,  the  method  of  ar¬ 
ray  averaging  and  differencing  are  significantly  simpler  than  the  bootstrap  based  method. 
The  computational  cost  of  computing  both  an  array  average  and  an  array  difference 
is  comparable  to  computing  the  test  statistic  (3).  However,  bootstrapping  requires  re¬ 
computation  of  this  statistic  many  times  and  thus  implies  a  larger  computational  burden. 
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VII.  conclusion 


Selecting  of  auto-  and  cross-terms  from  an  observed  TFD  has  been  cast  here  as  a  detec¬ 
tion  problem.  A  detection  approach  allows  performance  measurement  and  comparison  of 
various  schemes  via  the  operating  characteristic  curves.  The  idea  of  array  differencing  has 
been  introduced  and  combined  with  the  method  of  array  averaging  to  form  the  basis  of  a 
point  selection  algorithm.  Performance  of  this  scheme  has  been  evaluated  and  compared 
to  a  bootstrap  based  method,  showing  both  performance  improvement  and  computational 
savings. 
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WVD  of  Sensor  1 


Fig.  1.  (a)  WVD  of  chirp  signals  at  the  first  sensor,  (b)  Array  differenced  WVD. 
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Fig.  2.  Point  classification  scheme  for  selection  of  auto-  and  cross-terms  from  AA  or  AD  TFD’s 


OC's  for  auto-  and  cross-terms  detection 
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Fig.  3.  OC’s  of  methods  one  and  two. 
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Fig.  4.  OC’s  of  methods  one  and  two. 
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